_CR-1980*° 


j y - 07 ~c /s. 


NASA 



Supersonic Quiet-Tunnel Development (A I 7> 
for Laminar- Turbulent Transition Research: 

Final Report 

for NASA Langley Grant NAG- 1-1 607 


Steven P. Schneider 
Assistant Professor of Aerodynamics 
School of Aeronautical and Astronautical Engineering 

Purdue University 
West Lafayette, IN 47907-1282 



n ?) r o a 


ft m 

> C X 


O' CO 

3C i~t > 


• m 

Hm w 


> 

2? H > 


M TO 

> 1 I 


O 

PO H O 


O X 

1 C 7* 



H Z \ 


•n c z m 



73 m ^ 


— 3 

oo r oo 


T> ft 

c o 


C — 

r o ^ 


1 

m m o 


a z < w 


C ft 

h m 


ft T) 

r 


0 

H O CO 


c 

70 C 


3 f* 

> ac 


m, <m 

z m m 


< 

73 


* 

H H i/1 


^ ft 

H O 


CJ 

•-« -n z 


• 

o o ^ 





ru O 






1 


o 



UJ 



V. 



o 



NO 



o 

C 

z 

o 

3 


* 

O 

KJ\ 


— * 

\ 

UJ 

ft 

N) 



* 

* 


UJ 



o 



fNJ 


Abstract 

This grant supported research into quiet-flow supersonic wind-tunnels, be- 
tween February 1994 and February 1995. Quiet-flow nozzles operate with lam- 
inar nozzle-wall boundary layers, in order to provide low-disturbance flow for 
studies of laminar-turbulent transition under conditions comparable to flight. 
Major accomplishments include: (1) development of the Purdue Quiet-Flow 
Ludwieg Tube, (2) computational evaluation of the square nozzle concept for 
quiet-flow nozzles, and (3) measurement of the presence of early transition 
on the flat sidewalls of the NASA LaRC Mach 3.5 supersonic low-disturbance 
tunnel. Since items (1) and (2) are described in the final report for companion 
grant NAG-1-1133, only item (3) is described here. 


1 Introduction 

The general goals and background for this program were reported on recently 
in the final report for NASA Grant NAG-1-1133 [7], so these will not be 
repeated here. The grant reported on here supported additional graduate- 
student efforts closely related to those supported under NAG- 1-1 133. Items 
(1) and (2) from the abstract were reported on recently in the final report for 
NAG-1-1133 [7]. Some of the details for item (2) are contained in reference 
[1], which was completed after submission of reference [7] and is appended to 
this report. No further summary of items (1) and (2) will be presented here. 
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Grant NAG- 1-1607 also supported hot-wire measurements in the boundary 
layers on the flat sidewalls of the 2D NASA LaRC Mach 3.5 quiet tunnel, 
which were carried out during the summer of 1994 by Christine Haven (who 
was also supported under NAG-1-1133). These measurements were carried out 
in order to determine the Reynolds numbers at which the boundary layers on 
the sidewalls became turbulent. Computations of the boundary layers on the 
sidewalls of square nozzles [7] had provided clear evidence supporting earlier 
conjectures that transition on the sidewalls of square and rectangular nozzles 
is dominated by the crossflow instability. These computations also indicated 
that transition would occur very early due to this crossflow instability. This 
expectation of early transition on flat sidewalls subject to the 3D crossflow 
instability was supported by earlier unpublished pitot-probe measurements on 
the flat sidewalls of the NASA LaRC Mach 3.5 supersonic low-disturbance 
tunnel (SLDT) (S. Wilkinson, private communication). 

This early transition on the flat sidewalls is very significant to design of 
quiet-tunnel nozzles, for it suggests the following: 

1. The bleed slots for the flat sidewalls probably do not affect transition on 
the sidewalls, since transition still occurs near the throat due to crossflow. 
The turbulent boundary layer will be somewhat thinner with the bleed 
slots present, but it is hard to see why this would have a substantial 
effect on the size of the quiet flow region, although it would affect the 
radiated frequencies. Since these bleed slots add substantial mechanical 
complexity and cost, they should probably be omitted from future 2D 
designs. 

2. The rapid-expansion design is essential to the success of the Mach 3.5 
SLDT (with the 2D nozzle [2, 3]). The size of the quiet-flow region in 
a 2D nozzle is probably dominated by the noise radiating inward along 
Mach lines from the sidewalls (beginning at a point where the sidewall 
Mach number is 1.5 to 2, and the radiation becomes significant [5]). 
Inviscid computations carried out with the codes described in reference 
[7] indicate that longer nozzles cause the noise from the sidewalls to 
affect the centerline flow closer to the beginning of the uniform flow 
region. This makes a rapid expansion crucial to limiting the sidewall 
problem, since the only alternative is a very wide nozzle in which most 
of the massflow is wasted. However, recent work has found that a slow- 
expansion design is essential to achieving high Reynolds number quiet 
flow through reductions in the Gortler instability on the curved nozzle 
walls [9]. This conflict between the requirements for limiting the noise 
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from the sidewalls and the curved walls makes 2D nozzles unattractive 
for future high-Reynolds number quiet-flow designs. 

The significance of these sidewall transition issues made further measure- 
ments of transition on the sidewalls of the existing nozzle a priority. Clear 
measurements of the transition location would remove some of the uncertainty 
involved in the arguments. In order to obtain these measurements, Christine 
Haven spent the summer of 1994 in residence at NASA Langley attempting to 
set up and perform hot-wire measurements of the state of the sidewall bound- 
ary layer under various conditions. These measurements form the second part 
of her M.S. thesis [4]. Less than two days of tunnel access were provided, 
due to other priorities for the tunnel schedule. Unfortunately, this proved to 
be a severe limitation, since unexpected difficulties with noise generated by 
the traversing mechanism seriously corrupted the hot-wire data, and available 
resources did not allow for correcting the problem and re-acquiring the data. 

In spite of these difficulties with the noise in the data, some fairly definite 
conclusions can be drawn regarding the state of the sidewall boundary layer. 
The evidence for these conclusions is summarized in the balance of this report. 


2 Hot-Wire Measurements in the Flat Side- 
wall of the Mach 3.5 SLDT 

Most of the details of Christine Haven’s work in the SLDT at Langley are 
contained in her M.S. thesis [4]. Since some 256 megabytes of data were 
acquired, not all aspects of this work can be presented. The thesis has not 
been appended to this report, since the noisy quality of the data make the 
conclusions less clearcut, and the presentation does not match the quality of 
reference [1]. 

The data was acquired with a LeCroy model 9424 oscilloscope with 8-bit 
resolution. Records of 40,000 bytes were acquired at 400kHz using a prototype 
constant-voltage anemometer made by AS&M. The wire was placed about 
0.025 inches above the flat wall of the nozzle, and was moved streamwise and 
spanwise to sample different locations, which were mostly in the corners and 
in the exit plane. Data was also taken on the centerline of the tunnel for 
comparison to earlier measurements. Further details can be found in reference 

[4]- 

The large amount of noise in the data is evident from Figure 1, which 
shows the power spectrum of the hot-wire data at the nozzle exit. The 
conditions for the 3 files plotted can be found in Table 1. File rl_45 was 
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Figure 1: Power Spectrum of Hot-Wire Data at Various Conditions 


obtained under what were probably quiet flow conditions, whereas the other 
two files are probably of turbulent flow. The data show many spikes which were 
traced to electromagnetic interference from the DC motors in the traversing 
mechanism. Since most of the noise is above 50kHz, the signal to noise ratio 
can be enhanced by looking only at the signal content in the 0-50kHz band. 

Haven plots the power content in this 0-50kHz band in Table 2 in her thesis 
[4], and in figures 38-45. The power content is really in units of volts squared, 
but is mislabeled there as volts squared per second. In addition, there is an 
erroneous constant factor in her results. The original data that is cited in 
Table 2 of the thesis was reprocessed using a carefully checked spectral anal- 
ysis program written by the author and based on reference [6]. This spectral 
program computes and checks Parseval’s theorem for each dataset processed. 
The results are shown in Table 1. 

Here, x is the streamwise distance from the throat, in inches, y is the 
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y 

Pi 0 

Vm 

V 

v 50, h 

V'so 

Vso./i/Wo 

VsoAm 

dl.45 

3.95 

1.92 

11.1 

2.24 

2.35 

0.276 

8.5 

0.12 

jl-1 

15.07 

3.01 

34.0 

2.52 

4.10 

0.483 

8.5 

0.19 

kl.l 

15.47 

3.00 

4.8 

2.08 

0.73 

0.087 

8.4 

0.04 

pl-1 

15.47 

3.00 

28.2 

2.38 

1.54 

0.182 

8.5 

0.08 

rl.45 

4.22 

0.0 

9.1 

2.61 

0.09 

0.011 

8.2 

0.004 

xl.l 

15.22 

0.0 

68.3 

3.05 

0.30 

0.036 

8.3 

0.012 


Table 1: Sample Data from Haven (1995) 


spanwise distance from the centerplane, in inches, Pq is the stagnation pressure 
in psia, and V m is the mean CVA output voltage in volts. Files dl_45 and j 1_1 
were measured in the corner, kl_l and pl_l were measured on the sidewall in 
the exit plane, and rl.45 and xl_l were measured on the tunnel centerline. 
Also, K 5 ' 0 h is the rms power content in the fluctuating signal between 0 and 
50kHz, as computed by Haven (the square root of the data in column 6 of her 
table 2, p. 62 in [4]), and V£ 0 is the rms power content in the fluctuating signal 
between 0 and 50kHz as recomputed by the author. This makes V$ 0 /V m the 
normalized power content in the 0 to 50kHz range. Note that the ratio of V$ 0 k 
to V ' 0 is consistently about 8.4. It thus appears that Haven has an erroneous 
but consistent constant factor of 8.4 in her conversion from the signal to the 
integrated power. The rest of the data reported here will use Haven’s processed 
results, but correct them using this factor of 8.4. 

Since the signals and spectra are corrupted by noise, it is difficult to de- 
termine which of the data correspond to laminar conditions and which to 
turbulent. Without knowledge of the amplitude of the fluctuations, only a 
rather uncertain assessment can be made [4]. However, Haven acquired data 
on the centerline of the tunnel at known freestream conditions with the same 
hot wire used to acquire most of the rest of the data. The massflow can be 
computed for these conditions and compared to the hot-wire mean voltage. 
The results are presented in figure 2, which reproduces figure 19 from refer- 
ence [4]. It can be seen that the data form a consistent set. Since all the 
data was acquired with the same overheat ratio, 1.3, which is fairly high, it 
is reasonable to assume that the voltage is primarily responding to massflow 
fluctuations. Since the stagnation temperature was not changed during the 
run, and neither was the hot-wire voltage, the only mean temperature change 
would be that present in the boundary layer due to viscous heating effects. 
Since this cannot be determined from Haven’s data, which was obtained at a 
single overheat, it will be neglected here. The assumption will be made that 
the wire is responding primarily to massflow fluctuations, due to the high over- 
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Figure 2: Calibration of Wire One Used by Haven 

heat. This is a common assumption in CTA work [8], although the validity in 
the present case is unknown. 

If the voltage V is then assumed to depend only on the massflow pu, we 
can derive a relation between the normalized voltage fluctuations and the 
normalized massflow fluctuations by assuming small disturbances. Using a 
Taylor expansion, we find that 

v'/v = _ 

' m d(pu) V m ( pu) m 

Here, V‘ is the fluctuating part of the voltage, V m is the mean, (pu)' is the 
fluctuating part of the massflow, and ( pu) m is the mean. Although the small 
disturbance approximation is not valid for some of the data, it is a useful 
first approximation. Examination of figure 2 shows that dV/d(pu) varies from 
roughly 1 at V m = 2.6 to about 1/4 at V m = 3.0, and ( pu) m /V m varies from 
about 1/20 at V m = 2.6 to about 1/4 at = 3.0. The combined conversion 
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factor, [dV/d(pu)][(pu) m /V m ], varies from about 1/20 to about 1/16 over the 
same range. Values of this factor for smaller values of V m will be smaller than 
1/20, by an unknown amount. This makes the normalized massflow fluctua- 
tions about 16 to 20 times larger than the normalized voltage fluctuations, for 
small fluctuation levels. 

Using this rough calibration, we can then estimate the values of the normal- 
ized rms massflow fluctuations using the normalized rms voltage fluctuations 
shown in Table 1. It seems clear that only file rl_45 on the centerline can 
represent non-turbulent fluctuations. The fluctuations present even in this file 
are much larger than those measured in earlier, accurate LaRC data, but the 
difference is probably due to the additional electromagnetic interference noise 
in the present data. The data in the boundary layer from the flat sidewall all 
contain normalized massflow fluctuations that seem to be above 10 percent, 
and the corner data also seem to be above 10 percent. This amplitude data 
provide additional evidence that the flow on the flat sidewall of the nozzle is 
turbulent, although the several assumptions present in the calibration preclude 
drawing a firm conclusion. 

Data for the centerline fluctuations reduced using the process is shown 
in figure 3. These were taken under conditions where previous LaRC data 
showed a transition from quiet to noisy flow at about 11 inches from the throat. 
Since the mean voltage is about 3 volts, the normalized massflow fluctuations 
can be taken as about 16 times the normalized voltage fluctuations. The high 
fluctuation levels present in the upstream regions are presumably due to the 
electromagnetic noise from the traverse. 

The rms fluctuation levels for the lowest Reynolds number sidewall data 
are shown in figure 4. The data were taken at a total pressure of about 9 psia, 
the lowest pressure achieved during the runs. Four lines are plotted, for data 
taken along four different streamwise lines that are different distances from 
the curved wall. The four values of y ex i t for the four curves are the distance 
between the data-acquisition line and the centerplane, at the exit. The line 
with the smallest value of y ex i t is still fairly near the corner, about 6/7 of the 
distance from the centerplane to the corner. The figure reproduces the data in 
Figure 39 of reference [4], except for the corrected scaling. It appears from the 
large values of the normalized voltage fluctuations that the flow is turbulent 
for all cases. The data for the y ex it — 2.63 case near the throat exhibits a 
large amplitude spike in the spectra at about 15 kHz; this spike decreases 
in amplitude as the probe is moved downstream, until it is washed out in 
the general turbulence. Although this interesting spike may represent the 
instability waves that cause the transition, it is not emphasized here because 
there is no way to be sure that it was not caused by probe vibration or some 
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Figure 3: Normalized Fluctuations on the Centerline 
other interfering effect. 

3 Summary 

This graduate-education grant supported in part work by three master’s stu- 
dents: Timothy Alcenius, Christine Haven, and Laura Randall. Although 
Alcenius was interested in continuing for a PhD, the early termination of the 
original 3-year graduate-student-education grant (after one year) forced Alce- 
nius to obtain other employment. His work, and the first part of Haven’s work, 
is reported on elsewhere [7]. The second part of Haven’s work involved hot- 
wire measurements in the boundary layer of the flat sidewalls in the NASA 
Langley Mach 3.5 low-disturbance tunnel. As shown above, these provide fur- 
ther evidence that transition occurs very early on these flat sidewalls, which 
are in effect always turbulent. Seven months of Randall’s two-year M.S. pro- 
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Figure 4: Normalized Fluctuations on the Flat Sidewall 

9 psia total pressure, various distances from centerplane 

gram were also supported by this grant, at half of the usual student stipend. 
Her M.S. thesis is expected to appear in July 1995. 
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ABSTRACT 


Alcenius, Timothy John. M. S., Purdue University, December, 1994. Development of 
Square Nozzles for High-Speed, Low-Disturbance Wind Tunnels. Major Professor: 
Steven P. Schneider. 

Wind tunnels with low freestream disturbances are required to advance boundary- 
layer stability and transition research and provide accurate transition prediction methods 
for advanced aircraft. Developments in quiet-tunnel technology are reviewed. Also 
reviewed are the types of nozzles used for wind tunnel designs and the boundary-layer 
disturbances inherent in these designs. Three-dimensional nozzles may be best suited for 
quiet-flow tunnels if the boundary-layer crossflow remains small. To study the three- 
dimensional design for quiet nozzles, the flow fields in two Mach 2.4 nozzles were 
numerically simulated. Some design results for a Mach 8 nozzle are also presented. 
Preliminary estimates of crossflow-induced transition in the wall boundary-layers of the 
Mach 2.4 nozzles indicate that transition is likely to occur near the throat in both nozzles. 
A simple analysis of the crossflow pressure gradient in the throat is also presented. This 
preliminary analysis indicates that a radius of curvature of 10 throat radii or greater may 
be necessary to eliminate transition in the throat region. 
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CHAPTER I INTRODUCTION 


Accurate methods of estimating the location of boundary-layer transition are 
required for the development of aircraft for the Twenty-first Century, such as the High- 
Speed Civil Transport and NASP. To verify advanced computational fluid dynamics 
(CFD) codes used in the design of such aircraft, new ground test facilities must be 
developed. These facilities must be capable of obtaining transition Reynolds numbers 
on models comparable to those from flight data. 

The many factors involved in high-speed boundary-layer transition make it 
difficult to understand (see Morkovin [1]). To advance current knowledge, stability 
experiments in both the supersonic and hypersonic regimes are needed to clarify 
experimental and theoretical discrepancies [2,3]. New low-disturbance wind tunnels are 
essential for this purpose [4], 

Pate and Schueler [5] conclusively showed that wind tunnel noise was the major 
influencing factor in the scatter of measured transition Reynolds numbers on similar 
experiments in different experimental facilities. The high noise levels generally 
observed were caused by eddy Mach wave radiation from the turbulent boundary-layers 
on the nozzle walls [6]. Efforts by NASA Langley Research Center to develop quiet 
supersonic wind tunnels have shown that laminar boundary-layers are required on the 
nozzle walls. Measured transition Reynolds numbers on test models in quiet tunnels 
were then found to be as high as those observed in flight [7]. 

To introduce the laminar-flow control problem for wind tunnels, quiet-tunnel 
developments will first be reviewed. Then the types of nozzles used in quiet-tunnel 
designs will be discussed along with their inherent advantages and disadvantages. 
Disturbances that can cause boundary-layer transition will then be reviewed with 
emphasis on their application to quiet-tunnel design. Finally, the objectives of this 
research will be presented. 
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1.1 Quiet-Tunnels 

The considerable scatter in transition data from various conventional tunnels 
showed that new tunnels with low-freestream disturbances were necessary to provide 
accurate transition measurements. Since Laufer [6] had shown that wind tunnel noise 
was primarily due to the noise radiated along Mach lines from turbulent nozzle-wall 
boundary-layers, one of the main goals was to keep the boundary-layers laminar as far 
downstream of the throat as possible. Attempts to use shielding devices [8] and rapid- 
expansion nozzles with fully turbulent wall boundary-layers [9] were made. However, 
maintaining laminar boundary-layers along the nozzle walls proved to be more practical 
than either of these methods. Reviews of the evolution of high-speed quiet tunnel 
technology can be found in references 4 and 10. 

The first operational quiet tunnel was the Mach 3.5 Pilot Tunnel at NASA 
Langley. Measurements showed that the quiet-test length in this tunnel decreased from 
25 to 13 cm as the freestream unit Reynolds number (Roo) ncreased from 1 to 3 x 10? 
nr* . Efforts to maximize the quiet length Reynolds number (Rax) were then made to 

advance future designs. The first theoretical advance was determining how the nozzle 
design parameters affected the length of the quiet-test region. Calculations using the e N 
method [11] showed that transition was caused by the formation and amplification of 
G Order vortices along the concave part of the nozzle walls, instead of by the 
amplification of Tollmien-Schlichting waves. This result led to the design of a "slow- 
expansion" nozzle [12]. 

The maximum wall angle at the inflection point of slow-expansion nozzles is 
much smaller than for rapid-expansion nozzles [11], With the use of small wall angles, 
it is possible to insert a region of radial flow corresponding to a straight line wall 
section, inclined at the wall angle, upstream of the inflection point. This design change 
delayed the onset of the GOrtler vortices and reduced their growth rates by moving the 
inflection point far downstream where the boundary-layer is much thicker and the 
radius of curvature of the concave wall is much larger. Even though the Gdrtler number 
doesn't change much, longer nozzles give smaller values of integrated amplification for 
the Gdrtler instability [13]. The Rax values are then considerably larger than in the 
rapid-expansion nozzles. Measurements confirmed these predictions. 

Reviews of the current facilities at NASA Langley are given in references 10 
and 14. Background for new facilities at NASA Ames and Purdue University is given 
in references 15 and 16 respectively. Heating of the nozzle walls can also be used to 
maintain laminar boundary-layers. Background for the nozzle used in the investigation 
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of this method at Montana State University is given in reference 17 and the heating 
results in reference 18. 


1.2 Nozzle Types 

Generally, three types of nozzles are used for supersonic wind tunnels. These 
are: axisymmetric (circular cross-sections), two-dimensional (two contoured walls and 
two flat sidewalls), and three-dimensional (four contoured walls with square or 
rectangular cross-sections). Each has advantages and disadvantages to its use in a quiet- 
tunnel design. These are discussed below. 


1.2.1 Axisymmetric Nozzles 

The main advantage in using an axisymmetric nozzle is that stability issues that 
are inherent in two and three-dimensional nozzles are not present. Specifically, flow in 
comers and boundary-layer crossflow, which can cause stability problems, do not occur. 
Also, this type of nozzle will generally require a smaller mass flow rate than a two- 
dimensional nozzle with the same throat height. Finally, tolerances are easier to 
maintain in hypersonic nozzles since they have small, heated throats. 

There are, however, several disadvantages in using this type of design. First, 
machining imperfections can cause disturbances that are then focused along the axis of 
the nozzle. This may require smaller machining tolerances during fabrication to assure 
that freestream noise levels would be small. Also, the surface finish is difficult to 
maintain in the most critical location, the throat, due to access problems. Another 
problem is that long nozzles must be built in sections. This means that the joints 
between sections must be perfect or the same focusing effect caused by wall waviness 
would occur due to the surface discontinuities. The main disadvantage, however, is that 
circular walls complicate the use of optical grade windows. This means that standard 
flow visualization techniques such as Schlerien would be difficult to use in this type of 
nozzle. Since optical diagnostics are a highly desirable feature, this is a major limiting 
factor. However, new focusing Schlerien systems make straight walls non-essential for 
flow visualization using this technique. 
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1.2.2 Two-Dimensional Nozzles 

The two-dimensional design has several advantages over the axisymmetric 
design. First, since all the walls are fabricated separately and the cross-sections are 
rectangular, disturbances created by machining imperfections will not focus on the 
nozzle centerline. The disturbances will instead be laterally distributed across the flow. 
For this reason, machining tolerances are not as strict as in the axisymmetric design. 
Second, flat sidewalls allow for the use of optical diagnostics. Third, the nozzle can be 
taken apart to allow access to the throat region for polishing and maintenance. Finally, 
the long span between the two sidewalls can accommodate wider test models than the 
axisymmetric or three-dimensional designs would allow. However, noise from the 
sidewall boundary-layer could cause premature transition on the model if the model is 
too wide. 

The main problem with two-dimensional nozzles for quiet-flow applications is 
that transition on the sidewalls and in the comers is not well understood. To avoid these 
problems, the nozzles are made wide enough so that noise radiated from the sidewall 
boundary-layers and the comers will not affect the test region. The drawbacks to this 
approach are the increased cost necessary for a large nozzle, and the larger mass flow 
rates required. Finally, in high Mach number, two-dimensional nozzles the throat 
becomes so narrow and hot that it becomes impossible to maintain tolerances. 


1.2.3 Three-Dimensional Nozzles 

Three-dimensional nozzles seem to provide a compromise between the other 
two designs. The large mass flow rates that are required for wide two-dimensional 
nozzles are not needed for this design (assuming the throat heights are the same). Also, 
since all the walls are machined separately, the focusing effect that is evident in the 
axisymmetric nozzle does not occur. Like two-dimensional nozzles, these nozzles can 
also be taken apart for access to the throat for polishing and maintenance. Finally, use 
of optical grade windows is easier than in the axisymmetric design. 

The main drawback in this type of nozzle is the lack of understanding of the 
transition mechanics due to crossflow and comerflow. If further investigations can 
show these effects to be minimal or provide reasonable measures to limit the growth of 
these disturbances, this type of nozzle design may prove to be the most viable for quiet - 
flow applications. 
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1.3 Types of Disturbances 

Transition is a complicated process which is initiated by disturbances entering 
the boundary-layer. The process of entering the boundary-layer is known as receptivity 
and presently is poorly understood. After entering the boundary-layer, disturbances can 
grow or decay in either a linear or nonlinear fashion depending on their size and how 
they interact. Finally, they become large enough that the boundary-layer begins to 
become unstable and turbulent bursts are formed. Soon after, the entire boundary-layer 
will be turbulent. Reshotko [2] reviews the current knowledge of receptivity and 
disturbance evolution. 

Many different types of disturbances exist that can adversely affect boundary - 
layer stability. Since quiet-flow nozzles require laminar boundary -layers along the 
nozzle walls as far downstream of the throat as possible, it is desirable to keep all 
disturbances to a minimum. The boundary-layer stability on wind tunnel walls is 
evaluated by analyzing first and second mode viscous instability waves (first mode is 
otherwise known as Tollmien-Schlichting waves), Taylor-GOrtler vortices, wall 
roughness effects, comer flows, and crossflow. The effects that these disturbances have 
on boundary-layer stability can be estimated using Reynolds number correlations, linear 
stability analysis (e^), or the Parabolized Stability Equations (PSE) [19, 20]. However, 
these techniques are not capable of predicting transition when the disturbances become 
nonlinear. This happens when two of more disturbance waves interact and create a 
wave which is nonlinear. 


1.3.1 Tollmien-Schlichting Waves 

The continuity and Navier-Stokes relations can be manipulated to obtain the set 
of linearized disturbance equations. By assuming incompressible flow, no buoyancy 
effects, no curvature, and parallel flow for the basic flow and the disturbances, the most 
general form of a three-dimensional disturbance is found to be a traveling wave whose 
amplitude varies with y and moves along the wall at angle <|> with respect to the x-axis 
[21]. These are referred to as Tollmien-Schlichting waves and can be written as: 


(u,v,w,p) = [u(y), v(y),w(y),p(y)]exp{i<x(xcos<|> + zsin<j) - ct)} 


( 1 ) 
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where 

a Wave Number 

c Phase Speed 

co Frequency (= ac) 

These waves form one type of laminar-flow instability. 

To gain insight into the effect of these waves on nozzle-wall boundary -layer 
stability, an understanding of how pressure gradient affects wave growth must be 
developed. Nozzles with supersonic test sections have favorable pressure gradients if 
no shock waves are present. Wazzan, as shown in White [21], computed the neutral 
stability curves for the Falkner-Skan wedge-flow profiles. The results of this work 
showed that favorable pressure gradients strongly retard the amplification of Tollmien- 
Schlichting waves while adverse gradients promote them. 

The amplification of TS waves must still be considered, however, especially in 
slow-expansion nozzles. The e^ method from linear stability theory is the preferred 
method for computing amplification of these waves. The value of Nts at the location 
where Nq = 9 (the location where transition takes place due to the GOrtler instability) 
was negligible for the Mach 3.5 rapid-expansion nozzle but increased for the slow- 
expansion nozzles to 2.3 for the Mach 3.5 axisymmetric nozzle [12], to 3.6 for the 
Mach 6 axisymmetric nozzle [22], and to 4.5 for the Mach 2.4 axisymmetric nozzle 
[23]. Although, TS waves have not become large enough to cause transition in slow- 
expansion nozzle designs, this must be verified for each new design. 


1.3.2 Taylor-GOrtler Vortices 

Taylor-GOrtler vortices are counter-rotating vortices in viscous regions created 
by centrifugal instabilities. Rayleigh, as quoted in White [21], found that: "An inviscid 
rotating flow is unstable if the square of its circulation decreases outward." In other 
words, if the square of the product of the radius and the velocity decreases with 
increasing radius, the flow is unstable. Figure 1 [24] shows these vortices in the 
boundary-layer flow over a concave wall. 

As Beckwith, Malik, and Chen [11] noted, in supersonic nozzles, a uniform test 
region cannot be produced without a concave wall to cancel initial expansion waves. 
Therefore, the GOrtler instability will be present in any supersonic nozzle design. It also 
should be noted that a concave region must also exist in the subsonic portion of the 
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Figure 1 - Gdrtler Vortices within a Boundary-Layer Flow over a Concave Wall [24] 


nozzle to bring the flow from the settling chamber to the nozzle throat. However, any 
instabilities generated in this region can be sucked away by the use of bleed slots just 
downstream of the inflection point. 

Since it has been shown that boundary-layer transition on the walls of most 
quiet-tunnel nozzles is due to the GtJrtler instability rather than Tollmien-Schlichting 
waves [11], the main goal in new nozzle designs is to delay and limit the growth of the 
Gbrtler vortices as was discussed in the section on quiet- tunnels. A straight wall 
section, inclined at the maximum wall angle, is included in new designs for this 
purpose. 


1.3.3 Surface Roughness 

Experimental and computational data analyzing the effect of wall roughness on 
transition is presently lacking [14]. Therefore, it is difficult to determine exactly how 
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disturbances created by machining imperfections will affect boundary- layer stability. 
Beckwith, Malik, and Chen [11] have provided, however, a limited set of data for 
experiments in the NASA Langley Mach 3.5 Supersonic Low-Disturbance Pilot Tunnel 
with the tunnel walls "clean", "repolished", and "dirty". The "clean" condition 
corresponded to the nozzle being wiped with a lint-free cloth dampened with alcohol 
followed by vacuuming to remove atmospheric dust and lint. The "repolished" 
condition contained a significant reduction in the rms and peak-to- valley profilometer 
readings due to the polishing. Finally, the "dirty" condition occurred when there was a 
visible build-up of lint, dust, and other contaminants before a run. Conclusions from 
these runs showed that the repolishing operation resulted in a significant increase in 
laminar flow in the wall boundary-layer. Therefore, walls with small machining 
tolerances, and frequent cleaning and polishing seem necessary for quiet-tunnel 
operation. 

As a compromise between performance and machining capabilities in more 
recent quiet-flow tunnels, finish tolerances at the throat with maximum deviations in the 
range of 2 microinch rms have typically been specified. 


1.3.4 Comer Flow 

Comer flow is a complex three-dimensional flow formed by two intersecting 
surfaces. The three-dimensionality is enhanced by the presence of crossflow along the 
walls that transport mass, momentum, and energy into and out of the comer region [25] . 
It has been shown that two counter-rotating vortices centered about the comer bisector 
are also present in laminar supersonic flow through a duct with a mild adverse pressure 
gradient [26]. Clearly, the different factors that contribute to corner flow make 
determining the flow stability a difficult problem. 

The boundary-layer in a streamwise comer has been investigated. Recent 
examples of these studies are given in references [25-28]. The only known 
experimental study of transition in a comer flow with a favorable pressure gradient is 
that of Misu as reported by Beckwith [29], It was shown that reverse transition in the 
corner of a rectangular shaped contraction occurred at about the same streamwise 
location as on the walls. Following these results, it would seem reasonable to expect a 
similar result for transition in square nozzles since they also contain favorable pressure 
gradients. Although experiments by Zamir [30] were conducted close to zero 
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streamwise pressure gradient, he concluded that for Reynolds numbers greater than 10^, 
the comer boundary- layer is only stable with some favorable streamwise pressure 
gradient. Again, it is hoped that the favorable pressure gradient which exists in square 
nozzles will be helpful in preventing transition in the comers. 

Pfenninger and Syberg [31] have found that premature transition in the comer 
boundary-layer at the juncture of aircraft components can be prevented by suction. This 
result has been verified experimentally by Goldsmith (cited in reference 29). These 
results suggest that suction at the comers of the square nozzles could be used if 
required. 


1.3.5 Crossflow 

Various fluid flows occur in a way that creates a pressure gradient normal to the 
mean flow direction near a body. This pressure gradient will cause the development of 
a velocity component inside the boundary-layer that is perpendicular to the local 
inviscid flow and parallel to the wall. This is known as crossflow. Such applications 
consist of but are not limited to: flow over swept wings, flow on rotating disks and 
rotating axisymmetric bodies, flow in corners, and attachment- line flows. The 
crossflow profile is defined by a maximum somewhere in the middle of the boundary- 
layer and goes to zero at the boundary-layer edge and at the wall. A schematic of a 
three-dimensional boundary-layer profile along with the streamwise and crossflow 
components is shown in figure 2 [32], The crossflow profile contains an inflection 
point which is known to be dynamically unstable. This instability often leads to the 
formation of corotating vortices, which can be stationary or traveling crossflow 
disturbances. 

Two and three-dimensional nozzles exhibit a transverse pressure gradient due to 
the uneven expansion along the walls. This gradient will cause crossflow in the 
boundary-layers of these nozzles. Since crossflow is poorly understood for nozzle 
flows, a thorough investigation of nozzle wall boundary-layer stability in the presence 
of crossflow must be completed before three-dimensional nozzles will be proven useful 
for quiet-tunnel designs. 

The phenomena of crossflow and its effects on transition has been studied 
extensively, and comprehensive reviews are given in references [32-37], Since a vast 
amount of material exists on the subject, only material relevant to this research will be 
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Figure 2 - Velocity Components within a Three-Dimensional Boundary-Layer [32] 


covered. First, a discussion of the crossflow Reynolds number and its use in estimating 
transition is presented. Then some basic fundamentals of crossflow in subsonic and 
supersonic flow will be discussed. 


L3 .5 . 1 Crossflow Reynolds Number 

The earliest attempt to evaluate boundary-layer stability and transition in the 
presence of crossflow was through the development of the crossflow Reynolds number. 
The equation for this is shown below: 

R . WmAx6|Q (2) 

v e 


where (see figure 3) [38] 

WmAX Maximum crossflow velocity 

5 10 Maximum distance from wall where crossflow is 10% of W MAX 

ve Kinematic viscosity at edge of boundary-layer 
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Even though linear stability analysis is the technique recommended for determining 
transition location, there is still a desire to have parameters based purely on basic-state 
profiles which can be correlated to transition location. These parameters can then easily 
be incorporated into existing programs and would require much less time and 
computing effort than e^ methods to find the location of crossflow dominated 
transition. For this reason, crossflow Reynolds number techniques for estimating 
transition location are still used for preliminary design purposes. 

The introduction of the crossflow Reynolds number is credited to Owen and 
Randall [32 and 39]. They also suggested that transition occurs when the crossflow 
Reynolds number becomes on the order of 150 as observed in two-dimensional wing 
experiments at different sweep angles (Poll [37] discusses the details of this 
correlation). Further review by Poll shows that although initial experiments agreed with 
this correlation fairly well, experiments by Boltz, Kenyon, and Allen using a different 
wing section but the same sweep angles gave a transition Reynolds number from 200 to 
250. Finally, Poll states that experimental investigations of flow over a rotating disc 
found a transition crossflow Reynolds number of 680, significantly different than on 
swept wings. 

Various applications of the crossflow Reynolds number have been attempted to 
find a method for estimating the transition location which works for many different 
cases. Beasley, as cited by Amal [34, 36], defined a new crossflow Reynolds number 
based on an integral of the crossflow profile. This definition of the crossflow Reynolds 
number is shown in equation 3 below: 


1 Tie 

R 52 = — | w dy (3) 

v e 0 


where 

ve Kinematic viscosity at edge of boundary-layer 

Tie Distance to the boundary-layer edge along the wall normal 

w Crossflow velocity 

Using experimental data, it was suggested that transition would generally occur when 
this parameter reached 150. However, as shown before and noted by Amal [36] and 
Poll [37], the use of a single parameter to estimate transition cannot be good for all 
situations. For this reason, Amal and Coustols [34, 36] developed two criteria, Cl and 
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Figure 4 - Cl Crossflow Criterion [34] 


C2, which they felt could more accurately estimate transition. The first method. Cl, is 
based on Beasley's crossflow Reynolds number and the streamwise shape factor. A 
curve fit of experimental data was used to give a transition crossflow Reynolds number 
for various streamwise shape factors. This is shown in figure 4. Although some scatter 
of data still exists in this correlation, it appears to be much better than assuming that 
transition occurs at a constant value of the crossflow Reynolds number. 

The second method, C2, was developed in an attempt to reduce the scatter seen 
in the Cl criterion. This method is based on results from linear stability theory instead 
of a fully empirical development. Here, each disturbance angle is checked from 0° 
(fully streamwise profile) to 90° (fully crossflow profile). Then the critical Reynolds 
number can be found and checked against the streamwise shape factor and the 
freestream turbulence level to determine if transition has occurred. Although this 
method is based on linear stability theory, the authors admit that calculations to find the 
critical Reynolds number are long and costly, and not as accurate as e^. This deters its 
use for initial design estimations. 

Reed and Haynes [39] discovered that a single value for the crossflow Reynolds 
number does not correlate with transition in the supersonic regime either. Initially, they 
found that Chapman and Pate had concluded that the crossflow Reynolds number 
correlated well with the location of transition. Experiments by King on a yawed cone, 
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however, showed that the crossflow Reynolds number did not correlate with transition 
location. By completing calculations on rotating cones at various conditions, Reed and 
Haynes were then able to show that the crossflow Reynolds number at transition varied 
from 250 to 450, much different than the 150 observed in the subsonic regime. 

After showing the wide scatter in the crossflow Reynolds number for yawed and 
spinning cones, Reed and Haynes [39] presented a new crossflow Reynolds number, 
shown in equation 4, that included corrections for both compressibility and wall 
cooling. 


RcfRAH - HL (4) 

v e 

The equations for the compressibility correction, H, and the wall cooling, L, are given 
in equations 32-34 in chapter 2.2.3. Calculating this new crossflow Reynolds number 
from their numerical data on spinning cones, two points were established. First, their 
computed incompressible value fell among the compressible values. Second, a 
relationship seemed to exist between the new crossflow Reynolds number and the 
maximum crossflow velocity. Checking this new crossflow Reynolds number on the 
experimental data of King, two more important points were noted. First, it appeared 
they could estimate transition by dividing the new crossflow Reynolds number by the 
maximum crossflow velocity over the edge velocity (WMAXAJe) in percentages. This 
gave a constant value at transition of 44.0 for quiet flows and 33.7 for noisy flows. 
Secondly, they noted that the values of the new crossflow Reynolds number were in the 
neighborhood of 150 when WMAXAJe was on the order of 3% which was consistent 
with the results of Owen and Randall. 

Godil and Bertelrud [40] also noted that although the crossflow Reynolds 
number at transition is on the order of 200 in subsonic flows, the transition crossflow 
Reynolds number increased with Mach number as was shown by results from linear 
stability theory. Using these results, transition was then estimated from the standard 
crossflow Reynolds number with a correction for compressibility. Checking the results 
presented for both the corrected crossflow Reynolds number and linear stability theory 
on a biconvex wing at Mach 6, the new crossflow Reynolds seemed to estimate 
transition close to where it is predicted by linear stability theory (N = 9 for transition). 
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1.3.5.2 Crossflow in the Subsonic Regime 

Crossflow and its effects on transition have been studied extensively in subsonic 
flows, especially over wings, for laminar flow control design. On swept wings, a 
chordwise pressure gradient near the leading edge is the cause for the development of 
crossflow. Experiments by Dagenhart [41] on a 45° swept wing showed the existence 
of both stationary and traveling crossflow waves. Linear stability predicted the 
traveling wave frequencies accurately but stability results were about 30 percent larger 
than the observed stationary frequencies. Further work by Radeztsky, Reibert, and 
Saric [42] showed that linear theory does accurately predict mode shapes and expected 
wavelengths for stationary crossflow vortices. However, linear theory does not 
accurately predict the growth rates for low-amplitude waves. They report that this is a 
further extension of the experiments of Bippes and Dagenhart that showed the failure of 
linear theory for very large crossflow vortices. Finally, Dagenhart [38] concluded that 
crossflow disturbance amplification rates are proportional to the maximum crossflow 
velocity when all other factors are held constant. 

Analysis of flight measured data by linear stability theory [43] applied the 
theory of separation of transition mechanisms, the study of two or more disturbances 
which occur simultaneously but are assumed to act individually in causing the 
boundary-layer to become unstable, and found reasonable agreement with data from 
other experiments. Pfenninger [44] noted that this separation into independent parts is 
physically acceptable as long as strongly amplified crossflow and Tollmien-Schlichting 
waves do not occur simultaneously. In another attempt to study non-linear 
disturbances, Malik, Li, and Chang [45] completed a numerical study of the flow over a 
swept cylinder. The interaction studied was between traveling and stationary crossflow 
vortices. They conclude that when the disturbance amplitude reaches about 4 percent, 
the stationary and traveling modes begin to depart from their linear values. Another 
important conclusion from this work is that although the initial amplitude of the 
stationary vortex was much larger than the traveling modes, as the flow moves 
downstream, the traveling mode becomes the same order of magnitude as the stationary 
mode. 

Another type of flow that has been considered is flow over a rotating disk. In 
this problem, flow moves axially toward the disk and as the three-dimensional 
boundary-layer builds on the surface, the fluid is cast off the edge like a centrifugal 
pump. This causes a stationary crossflow instability on the disk which spirals outward. 
A compilation of data has shown that traveling waves have a higher amplification rate 
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than stationary waves according to theory [32], The importance of secondary 
instabilities still remains unclear for this problem (see reference 33 for a discussion of 
the effects of secondary instabilities in crossflow induced transition). 


1-3.5. 3 Crossflow in the Supersonic Regime 

Few numerical or experimental results have been obtained for crossflow in the 
compressible regime. Although Wang, Herbert, and Stucker [46] have completed some 
Parabolized Stability Equations calculations for a compressible swept wing flow, they 
"...have not pursued this supersonic case because data for comparison are lacking." 
Dagenhart [38] suggests that it seems that incompressible flow theory as it relates to 
crossflow is a physically reasonable approximation although it is somewhat 
conservative. Finally, Amal [34] states that the crossflow instability is dominated by 
the inflection point properties which are not affected by compressibility very much. For 
this reason, he also believes that incompressible criteria can be used to estimate 
transition in the compressible region. 

Numerical investigations of crossflow over various bodies at supersonic and 
hypersonic speeds showed, however, that transition in these flows does not agree with 
the predicted location using the crossflow Reynolds number. Reed and Haynes [39] 
show that the standard crossflow Reynolds number has a spread on the order of 200% 
for their calculations of flow over a rotating cone. They also show that the experimental 
results of King for flow over a yawed cone give even larger spreads than they observed 
for spinning cones. Finally, Godil and Bertelrud [40] stated that the transition crossflow 
Reynolds number increases with Mach number which was verified by the use of linear 
stability analysis. 

Supersonic flow inside a square or rectangular duct has been investigated by 
various researchers. Numerical results from Davis, Gessner, and Kerlick [26] of 
laminar flow through a square duct with a mild adverse pressure gradient show the 
development of two counter-rotating vortices around the comer bisector. They Fmd that 
this secondary flow causes a bulging of the total pressure contours in the comer region, 
a characteristic that is distinctly different from crossflow investigated in an unbounded 
comer problem. Numerical investigations of compressible turbulent flow through a 
square duct [25] showed similar total pressure contours and crossflow profiles to those 
measured in incompressible square duct flow at approximately the same streamwise 
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location. A numerical investigation of turbulent flow in square Mach 4.7 and Mach 6 
nozzles [47] showed some interesting flow phenomena that could also be expected in 
slow-expansion nozzles. Both nozzles showed the development of counter-rotating 
vortices in the comers that were then shed to the centerline. At the exit, these vortices 
were greater than 30 percent of the spanwise length of the wall. However, it is stated 
that agreement is not very good compared with experimental data close to the wall. 
Mach number contours at the exit also show small structures near the corner that are not 
discussed in the paper. This is probably the same type of structures shown by other 
investigators in comer regions. 

Comparisons of predictions and experiments for a Mach 3 two-dimensional 
nozzle [17] showed the development of crossflow along the nozzle walls. However, in 
this investigation, the crossflow originated at the center of the contoured wall and 
moved toward the comer and then down the sidewall. This created one vortex in the 
comer instead of the counter-rotating vortices observed in square ducts or unbounded 
comer problems. Another important conclusion in this work is that the calculated 
Gdrtler instability was insufficient to cause transition, and the observed crossflow could 
have contributed to the transition process. No details were provided for the affect of 
crossflow on transition, however. 

Experiments have also been performed on swept cylinders at Mach 10 [48]. The 
results showed that although no transition occurred, streaks generated by the crossflow 
instability were seen. Linear stability analysis applied to the flow provided two 
significant conclusions. First, increasing the wall temperature seemed to have a 
destabilizing effect on the stationary crossflow disturbance. Second, the maximum N 
factor calculated was less than 2. This shows that no transition should have occurred, 
verifying experimental results. 

Clearly, more work needs to be done in this area for a complete understanding 
of crossflow- induced instability and transition in high-speed applications. 


1.4 Objectives 

The objective of this research is to numerically study crossflow in two Mach 2.4 
square nozzle designs to estimate the boundary-layer stability. Specifically, a short and 
a long Mach 2.4 nozzle will be studied to determine if length has an effect on crossflow 
stability. The initial design for a Mach 8 nozzle is also presented as an initial step 
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towards examining the effects of Mach number on crossflow. To get a general idea of 
crossflow instability and transition, three crossflow Reynolds numbers will be 
calculated and simple transition estimators using these parameters will be studied. 
These are: 1) The standard crossflow Reynolds number with the transition estimator 
given by Godil and Bertelrud [40], 2) The new crossflow Reynolds number presented 
by Reed and Haynes [39], 3) The crossflow Reynolds number of Beasley [34, 36] with 
transition estimated using the Cl criteria of Amal and Coustols [34, 36]. Although the 
crossflow Reynolds number is a technique with limited accuracy, alternate techniques 
such as e^ or the Parabolized Stability Equations require far more resources than are 
available at this initial point in the design phase. It should also be noted that no data 
exists, to the authors knowledge, for crossflow- induced transition on nozzle wall 
boundary-layers with which these more accurate techniques can be evaluated. 
Therefore, the crossflow Reynolds numbers are used as an initial attempt to discover 
how crossflow will affect the boundary-layer. 

The nozzles for this study are designed using previously developed design 
techniques that have been found to limit certain disturbance amplifications. As 
discussed earlier, slow-expansion nozzles have been shown to limit the growth of 
Taylor-Gdrtler vortices while not significantly enhancing Tollmien-Schlichting waves. 
Therefore, this design will follow previous work and use a large radius of curvature at 
the throat, a small expansion angle, and a long radial flow region. It will then be 
assumed that neither of these disturbances becomes large enough to cause transition. 
Further investigation of these disturbances must be completed in the future for a full 
understanding of the overall boundary-layer stability. 

The subsonic portion of the nozzle must also be dealt with to insure a laminar 
boundary-layer upstream of the throat. To accomplish this, a bleed slot will be inserted 
in the nozzle downstream of the subsonic wall inflection point [8]. The purpose of this 
slot is to remove the contraction-wall boundary-layer and create a new laminar layer 
starting on the slot lip. The slot will not be considered in this portion of the research, so 
the upstream boundary in these calculations will be taken where the slot lip would begin 
for each nozzle. 

Although no design procedures exist for comer flows, this issue still must be 
addressed for a complete picture of the problem. It has been shown that transition in the 
comer occurs at roughly the same streamwise location as on the walls [29]. Other 
results su gg est that suction could be used to limit transition in the comer regions if 
necessary. Therefore, stability of the comer flow will not be considered here. 
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CHAPTER H NUMERICAL SIMULATION TECHNIQUE 


Various computer programs for the design of the nozzles, the numerical 
simulation, and the data post-processing were either implemented or developed. A 
discussion of these programs and the design procedures used follows. 


2.1 Nozzle Development 

Since most wind tunnel nozzles in use are two-dimensional or axisymmetric, 
few design methods exist for three-dimensional nozzles. Beckwith, Ridyard, and 
Cromer [49] developed a method for designing a three-dimensional nozzle using an 
axisymmetric flow field. In this procedure, an axisymmetric design is calculated by the 
method of characteristics in the supersonic region. A cross-sectional shape is then 
selected at the exit and the streamlines passing through the exit shape are tracked 
upstream. This provides a three-dimensional nozzle with an arbitrary cross-section. 

In quiet-flow nozzles, the three-dimensional contour is required starting from 
the bleed-slot lip. Therefore, the contour must also be found in the transonic region of 
the nozzle. Using a series solution developed by Hopkins and Hill [50], an 
axisymmetric wall contour can be developed for the transonic region of the nozzle. The 
streamline traces from the supersonic portion of the nozzle can then be continued into 
the transonic region, completing the nozzle contour. 

If this method is applied to a nozzle with a square test section, the cross-sections 
will only be square at the exit and at the throat if the flow is parallel and uniform. All 
other cross-sections will "bowed" in or out between the comers. This will be discussed 
in further detail in subsequent sections. A correction can be made by assuming a 
constant area and then each cross-section can be "squared-off'. This approximation has 
been used for the 31-Inch Mach 10 Wind Tunnel and two lower Mach number scramjet 
testing nozzles at NASA Langley as reported by Beckwith [29], A boundary-layer 
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correction is also included in the design of these nozzles. This consists of the addition 
of the local boundary-layer displacement thickness to the inviscid contour. The 
boundary-layer thickness is computed using a two-dimensional method with the 
centerplane pressure distribution. Flow surveys in the Mach 10 nozzle have shown that 
the flow uniformity is excellent. This indicates that the assumptions used in the design 
were satisfactory. 

The subsequent sections discuss the application of this method for the design of 
both Mach 2.4 nozzles and the preliminary Mach 8 design. 


2.1.1 Wall Contour Design 

The inviscid design for the wall contour is performed using the method of 
characteristics technique developed by Sivells [51] in the supersonic region, and the 
method of Hopkins and Hill [50] in the transonic region. A brief description of the 
codes used in this stage of the design follows. 

The beginning of any nozzle designed using this method is in the transonic 
region. The method developed by Hopkins and Hill uses an inverse solution to find the 
body geometry given a velocity distribution along the axis. Assuming steady, 
irrotational, adiabatic, shock-free flow of a perfect gas with constant specific heats, two 
differential equations that satisfy continuity and irrotationality were found using 
asymptotic methods. The equations were then transformed to use the velocity potential 
and the stream function as the independent variables. The dependent variables were 
then expressed in terms of a power series with the accuracy dependent on its order. 
These equations can then be solved point-by-point as long as the Mach number 
distribution is specified along the nozzle axis. A code was written [52] to solve these 
equations developed by Hopkins and Hill. This code is described below. 

First, the streamline that corresponded to the wall contour needed to be 
determined. Since the equations were normalized with respect to the throat radius, the 
streamline that had a minimum at y =1.0 corresponded to the wall contour. This 
streamline was found by using an iteration process, over Hopkins and Hill's variable £ 
to find the location where dy/dx = 0 when y = 1.0. 

The Mach number distribution along the axis could then be specified as long as 
the initial and final Mach numbers were known. These values were initially arbitrary 
and therefore a method was developed to find than. The starting location was taken 
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where the angle of the wall contour streamline was the same as the approach angle, 0 a . 
This location was chosen since it corresponds to the inflection point of the wall in the 
subsonic region. The location is, therefore, assumed to be the preliminary location of 
the bleed slot in these designs. The starting location was found by solving for the flow 
angle, as the dependent variable, using a bisection method, until the angle 0 a was found. 
The corresponding Mach number on the reference line was then used as the initial Mach 
number for the calculation. The final location was taken at the throat, where x = 0, so 
that the code could easily be patched with the method of characteristics code. This 
location was also found by solving for the dependent variable x using a bisection 
iteration method. After this location was found, the corresponding Mach number on the 
reference line was set as the final Mach number. 

The supersonic region was started at the end of the transonic flow solution, 
slightly downstream of the sonic line. Starting from this point, a characteristic net is 
developed up to the characteristic that defines the beginning of the radial flow region. 
The radial flow region is then patched to the solution with the requirement that the first 
and second derivatives of the Mach number on the centerline be continuous. Finally, 
the rest of the nozzle was calculated using the method of characteristics beginning at the 
characteristic that marks the end of the radial flow region. A fourth order polynomial is 
fit between Mb and Me to assure a smooth wall contour. The coordinates of the wall 
were determined by integrating the mass flow rate along the characteristic lines to 
satisfy the conservation of mass principle. 

After the solution to the flowfield was found, a square cross- section is defined in 
the exit plane and the corresponding streamlines are tracked upstream. Since the entire 
flow solution is known in the supersonic region, each of the desired streamlines could 
be found by its corresponding inflection angle. The limits of the streamlines were 0w 
and 0 W / V2.0 corresponding to the comerplane and centerplane of the square contour. 
However, in the transonic region, each of the desired streamlines must be found as if it 
were considered to be the wall. Therefore, the streamlines from where y = 1 / V2.0 
(corresponding to the square nozzle wall in the centerplane [see figure 5]) to when y = 
1.0 must be determined. This gives all the streamlines for the square nozzle contour. 
The streamlines from the Hopkins and Hill formulation in the transonic region are 
matched to the streamlines produced by Swell's code, at the streamwise location of the 
throat (x=0). It is important to note that although Swell's code uses a different 
asymptotic formulation for the transonic region, the two formulations match smoothly 
at the throat where they are patched together. This is not surprising since both 
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formulations are asymptotic solutions valid near the throat. 

The first phase in developing the desired nozzles was to obtain an axisymmetric 
contour. The data required for this were the radius of curvature at the throat (Ra), the 
approach angle (0 a ), the expansion angle (0w), the Mach number along the centerline at 
the end of the radial flow region (Mb), and the exit Mach number (Me). To keep the 
Gdrtler instability as small as possible, a large radius of curvature, small expansion 
angle, and long radial flow region are necessary. These are reflected in the design 
parameters given in table 1. Note the radius of curvature at the throat is non- 
dimensionalized with respect to the throat radius, so this value is dimensionless. 

Nozzle 1 is the first axisymmetric design and is referred to as the short Mach 2.4 
design from here on. The second Mach 2.4 nozzle (nozzle 2 in table 1, and referred to 
as the long Mach 2.4 nozzle) was designed using a smaller expansion angle, a smaller 
approach angle, and a larger radius of curvature. The design of the Mach 8 nozzle 
(nozzle 3 in table 1) was selected to match the size of the Mach 8 Variable Density 
Tunnel (M8VDT) [10] as closely as possible. 


Table #1 -Nozzle Design Parameters 


Nozzle 

Ra 

0a 

0w 

Mb 

Me 

Po 

(psia) 

To 

(°R) 

Roo 

(ft- 1 ) 

mm 

1 

5.5 

42.42° 

5.44° 

2.39 

n 

100 

540 

3.182 
x 10 7 

3.65 

2 

6.22 

35.36° 

2.83° 

2.39 


100 

540 

3.182 
x 10 7 

3.65 

3 

3.3 

30.0° 

6.50° 

7.75 

8.0 

1000 

1380 

9.764 
x 10 7 

0.525 


Note that the parameters R a , 0 a , and 0w arc along the diagonal of the square nozzle 
which corresponds to the radius of the axisymmetric design (see figure 5). To find the 
parameters along the centerplane for any wall of the square nozzle, those given above 
must be divided by the square root of two. 
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The parameters Rg, 0a* 0w> Mg, and Me we entered into the program, 
described above, that uses Sivells implementation of the method of characteristics for 
axisymmetric inviscid nozzle design with the method of Hopkins and Hill for the 
transonic portion. The streamlines corresponding to a squwe in the exit plane were 
written to a file for post-processing to develop the final three-dimensional contour. 

Once all the streamlines required from the axisymmetric contour were known, 
the final squwe contour can be completed. Full descriptions of the programs used can 
be found in references 52 and 53. A program named SNODEC [53] was developed by 
Michael Moen as a post-processing code for the development of square nozzles from 
axisymmetric designs. This code reads in the streamline data and creates the final 
square contour. While reading the data, the code asks for the number of planes 
downstream of the throat that are desired. The documentation suggests that no more 
than 200 planes be taken so as not to exceed the limits of the interpolation routine. 
However, the number of planes chosen should be the same order of magnitude as the 
number of grid points in the streamwise direction to maintain accuracy in the grid 
generation program. Therefore, 200 planes we chosen for all cases examined. Since 
the finest grid examined only has 305 points in the streamwise direction, 200 planes 
should be adequate for the interpolation routine. After the streamlines we read in and 
the data is interpolated to uniform stations, the squwe contour can be generated. Some 
bowing will occur in each cross-section except at the exit and at the throat, if the flow is 
pwallel, during this procedure. Figure 6 shows the difference in the z-coordinate 
between the corner and center through the short Mach 2.4 nozzle. If the difference in 
the z-coordinate between the corner and the center were zero, the cross-section would 
be a squwe with a height of z. However, some bowing can be seen both inward and 
outwwd. This is caused by the streamlines not being pwallel at every location. This 
bowing creates a contour that is not completely squwe. To squwe the cross-sections, a 
numerical integration is completed at every streamwise location to find the wea and 
then the squwe root of the wea is taken to obtain the squwe height. 

The nozzle contour is then completed by the calculation and addition of a 
boundary-layer correction. This correction is made by the addition of the displacement 
thickness calculated assuming a two-dimensional flow along the nozzle centerplane. 
Specifying the pressure distribution, calculated in the SNODEC code, and the wall 
coordinates, a program written by Harris and Blanchwd [54] for solving two- 
dimensional and axisymmetric boundary- layers was used. Also necessary we the throat 
radius, Prandtl number, stagnation pressure, stagnation temperature, and universal gas 
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Figure 6-Wall Bowing Due to Streamline Tracking of the Square Geometry 


constant. Given this data, the displacement thickness is calculated and then added to 
the inviscid contour by use of the SNODEC code. The SNODEC code is also used to 
calculate the pressure gradient between the corner and centerplane which will be 
compared to the results from the Navier-Stokes solutions. The final nozzle contour is 
given in terms of feet. 

The techniques described were then used to create the square contours for the 
two Mach 2.4 nozzles and the Mach 8 nozzle. The stagnation pressure, stagnation 
temperature, and throat radius are also shown in table 1 for each nozzle. To account for 
the location of the bleed-slot in these nozzles, they are cut-off at an approximate Mach 
number of 0.5 as an initial estimation. The final three-dimensional outline for the short 
Mach 2.4 nozzle is shown in figure 7, the Mach 2.4 long nozzle in figure 8, and the 
Mach 8 nozzle in figure 9. All of the units have been converted to inches in these 
figures. As was stated, the Mach 8 nozzle was designed to match the operating 
parameters of the M8VDT as closely as possible. The overall length of the nozzle using 
these design parameters is approximately 1.75" shorter than the M8VDT and the 
distance along the diagonal at the exit is 0.047" larger. Although a grid was generated 
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for the Mach 8 nozzle, simulations were not completed on this nozzle. Further 
investigations will be completed on this design at a later time. 


2. 1 .2 Grid Generation 

The grid used in all nozzles is a simple body fitted grid with constant spacing in 
the streamwise direction and stretching in the normal directions. Two programs were 
written to generate the grids. They are listed in appendix A. 

The first of these two programs takes the final square contour from the *.csd file 
generated by the SNODEC code. The structure for the *.csd file is shown in appendix 
A. The purpose of this program is to read in the contour from this file, convert it to 
metric units, and then write the coordinates in a generalized form. This form is: the x 
location, the minimum and maximum y values, and the minimum and maximum z 
values. Also written at the top of the data file are the minimum and maximum x values, 
to set the grid spacing in the next program. One important note, a number 1 must be 
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added to the bottom of the *.csd file to signify to this program that the data set is 
complete. 

The second program is used to generate the grid for the Navier-Stokes 
calculation. First, the data set that was created by the first program is read. Then the 
user is prompted for the stretching constant which controls the stretching near the wall. 
The choices for the stretching variable will be described later. Finally, the program asks 
the user if Tecplot [55] data is desired. The data generated for use in Tecplot consists of 
three files. The first file contains cross-sectional data at various constant x locations, 
the second contains a plane where y is the maximum value, also a wall contour, and the 
third contains a plane where z is zero, also a centerplane. These can be used to evaluate 
the grid without the unnecessary work of converting the full binary file. 

The grid dimensions are then set in the IMAX, JMAX, and KMAX variables. 
Next are constants defined for the grid stretching. These variables are set in the 
program since they remain constant for most cases. The variable XBAR can be set 
between zero and one and determines the weight of stretching near the wall versus the 
amount of stretching near the centerline. For these grids, this variable is chosen to be 
zero so all the stretching is near the wall. The variable EX is an exponential that limits 
how the centerline stretching affects points near the wall and how the wall stretching 
affects points near the centerline. However, when XBAR is zero, this parameter plays 
no role in the calculation. Variable SFACTC is the stretching variable for the grid near 
the centerline. This also plays no role in the calculation when XBAR is zero. This 
stretching function was originally developed by Dr. John Korte at NASA Langley and 
has been adapted for use in this code. 

An interpolation routine is set up so that no matter how many points are in the 
input file, an even spacing of points in the streamwise direction is developed. For this, 
a linear interpolation is used between points in the input file if necessary. Therefore, 
the number of points used to generate the *.csd file in SNODEC should be 
approximately the same number of points used in the grid in the streamwise direction. 
Finally, all of the points are calculated and written to a binary file that can be read by 
the Navier-Stokes solver. A listing of both of these codes as well as the file generated 
by the first code for the Mach 2.4 nozzle is given in appendix A. 

Since the flow is symmetric about any centerplane in the nozzle, the solution 
only needs to be calculated in one quadrant. Four calculations have been completed on 
the two Mach 2.4 nozzles. Three were done on the short Mach 2.4 nozzle as a grid 
resolution study. The grids used were 137x41x41, 201x65x65, and 305x97x97 
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corresponding to 2/3 and 3/2 of the medium grid. These grids were chosen to be 
multigridable up to at least four levels. The stretching constant was defined for these 
grids to be 1.00005. This was found by developing numerous grids using the medium 
resolution case to find a stretching parameter that would place approximately 20 points 
in the boundary-layer at the exit. The boundary-layer thickness was taken from the 
results from the Harris and Blanchard boundary-layer code. Then the same stretching 
was used on the coarse grid. The stretching was also attempted on the fme grid, but the 
high aspect ratios generated using this stretching caused too much artificial dissipation 
in the Navier-Stokes solver and the code would not converge. To correct this, a slightly 
larger stretching factor of 1.0005 was used for the fme grid. Further studies of this grid 
showed that the ratio of points in the boundary-layer was still approximately 3:2 
between the high resolution and medium resolution cases. 

The results from the grid resolution study have shown that the medium 
resolution grid is adequate for the crossflow Reynolds number techniques used. For 
this reason, the medium grid is also used on the Mach 2.4 long nozzle. Since the exit 
dimensions are the same, the medium resolution grid with a stretching factor of 1 .00005 
was adequate to resolve the boundary-layer properties. The problem with using this 
grid is since the long nozzle is almost twice as long as the short nozzle, the spacing in 
the streamwise direction becomes almost twice as large in the long nozzle. However, it 
will be shown that this should be adequate to resolve the flow. 


2.2 Numerical Solution and Crossflow Calculation 
The purpose of this section is to describe the techniques used by the three- 
dimensional Navier-Stokes solver and the data post-processing program. 
Documentation for the Navier-Stokes solver, LARCK (Langley Algorithm for Research 
in Chemical Kinetics) [56], is not yet available since the code is still under development 
at NASA Langley. Only the portions necessary for this problem are described. 


2.2.1 Governing Equations 

The equations to be solved for all of the nozzles considered are the full three- 
dimensional compressible Navier-Stokes equations. These are written in non- 


30 


dimensional form below 


au 3 e aF dc 

dt dx 3y dz 
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where U, E, F, and G are the vectors 
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which are non-dimensionalized using the relations 
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The components of the shear-stress tensor and the heat-flux vector in non-dimensional 
form are given by 
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and the perfect gas equations of state are 
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P = (y - l)pe 
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w 
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and ReL is the Reynolds number based on a reference length is 
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To solve these equations, constant specific heats are assumed as well as a constant 
Prandtl number. The viscosity is found by using Sutherland's Law. The variables 
included in the above equations are as follows 


x, y, z 
u, v, w 
P 

P 

T 

y 

Et 

e 

Uoo 

a 

R 

P 

Qx, qy, qz 

T xx> T yy. "tzz 
T xy. T xz. "tyz 
Pr 

ReL 

00 


Principal directions in the Cartesian coordinate system 

Velocities in the x, y, and z directions 

Density 

Pressure 

Temperature 

Ratio of Specific Heats 

Total Energy 

Internal Energy 

Freestream or Reference Velocity 
Speed of Sound 
Universal Gas Constant 
Viscosity 

Heat-Flux in the x, y, and z directions 
Shear-Stress Tensor along Specific Faces 

Prandtl Number 

Reynolds Number Based on a Reference Length L 
Subscript Referring to a Freestream or Reference Condition 


To solve these equations, various parameters need to be specified. In nozzle 
flow cases, where the freestream conditions are stagnation conditions, another point 
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needs to be chosen for the reference location. The point chosen for these calculations is 
where the Mach number is one, or the nozzle throat in a one-dimensional analysis. The 
input to the program is the stagnation pressure, stagnation temperature, and stagnation 
density which can be converted to values at the throat using isentropic relations. Also 
necessary are the Reynolds number at the reference location, Prandtl number, reference 
length, reference Mach number (1.0 for these cases), and the freestream velocity 
components (0.0 for both cases since freestream is stagnation). The value used for the 
reference length is 1 .0 so this variable drops out of the non-dimensionalization of the 
Navier-Stokes equations. The input geometry is also non-dimensionalized by this value 
so the input coordinates need not be changed. The inputs used for both Mach 2.4 
nozzles are given below: 


Stagnation Pressure 

689480 Pa 

Stagnation Temperature 

300 °K 

Stagnation Density 

8.0079 kg/m 3 

Reference Mach Number 

1.0 

Reference Length 

1.0 m 

Reynolds Number 

1.0439 x 10 8 

Prandtl Number 

0.72 

Ratio of Specific Heats 

1.4 


Given these parameters, the Navier-Stokes equations could then be solved for 
the dependent variables p, pu, pv, pw, and Et. This is done using a cell centered, 
central differencing scheme in space and a five-stage Runge-Kutta algorithm to advance 
the scheme in time. Techniques to accelerate convergence to steady-state are also 
included. These are implicit residual smoothing and multigrid acceleration. Each of 
these will be discussed in the next section. Two assumptions are also made for these 
computations. First, laminar flow is assumed for the entire flow. Since the flow should 
be laminar at the bleed slot lip, corresponding to the initial location of the calculation, 
this is the correct assumption until transition, which is the location that is being sought. 
To save calculation time and computer storage, it is also assumed that all the viscous 
terms containing derivatives in the stream wise direction are neglected. This is known 
as the "thin-layer" approximation to the Navier-Stokes equations. Since the viscous 
terms should be relatively small in the streamwise direction and the grid spacing is 
relatively large in this direction, dropping these terms should not affect the final 
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solution. 


2.2.2 Numerical Algorithm 

The space discretization begins by assuming that all the dependent variables are 
known at the center of all the grid cells. The flux vectors are then evaluated by finding 
the average values at each cell face. Finally, the total flux through any individual cell 
can be evaluated. This technique is described in more detail in references [57] and [58]. 
These schemes reduce to second order central difference schemes on smooth Cartesian 
meshes. 

The equations are advanced in time by the use of a five-stage Runge-Kutta 
scheme. This approach is shown completely in reference [59]. The dependent variables 
at the (q+1) stage is given by 


W 


(q+i) = w (o ) _ ttq+1 ^ Q |w(qyj _ ^ YqrD (w( r ' 
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= Q c (V q ^ +Q d (w(°)) 


where 


w(°) 

Solution at time level n 

aq+1 

Tunge-Kutta Coefficients 

At 

Time Step 


Weighing Factors for the Artificial Dissipation 

a 

Cell Volume 

Q c (w(q)) 

Convective Fluxes at time step q 

Qd(W(0)) 

Dissipative Fluxes at time level n 

D(W( r )) 

Artificial Dissipation 


For a five-stage scheme, the variable q varies from zero to four. When q+1 is five, this 
corresponds to the new time level n+1. The coefficient alpha has been determined such 
that the scheme has the largest hyperbolic stability limit. The coefficients used for these 


cases are 
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a i = 1/4 a2 = 1/6 a3 = 3/8 04 = 1/2 as = 1 

This scheme exhibits the good high-frequency damping that is required for a rapidly 
convergent multigrid method [59], It should be noted that the dissipative fluxes are 
only evaluated on the first stage and appear not to have an affect on the scheme stability 
[59], This allows for a significant reduction in computational effort for the scheme. 
For a good parabolic stability, the artificial dissipation terms are evaluated on the first, 
third, and Fifth stages. The weighing factors for the dissipation terms must satisfy the 
condition 


lYqr=l (12) 

The individual factors are defined in reference [59]. The maximum CFL number 
allowed by stability for this scheme is 3.75, which is about 30 percent larger than for 
the fourth-order Runge-Kutta scheme with one evaluation of the artificial dissipation. 
Finally, the steady-state solution is independent of the time step, and therefore, 
convergence acceleration techniques are easily employed in this method. 

Three methods are used to accelerate convergence of the basic scheme. These 
consist of local time stepping, implicit residual smoothing, and multigrid. The smallest 
time step computed in a computational domain is used for all cells in a global time 
stepping method to advance the calculation in time. However, this method is inefficient 
in the cells where a larger time step would be allowed by stability. Local time stepping 
is therefore applied to march the solution to steady-state using the maximum time step 
allowed by stability in each individual cell. Using this method, steady-state will be 
reached in some cells more quickly than by using a global time stepping method. 
Therefore, steady-state solutions will be available for flux calculations in cells with 
smaller stability limits, allowing for faster overall convergence. Both convection and 
diffusion stability limits are included in the calculation of the time step in any cell. 

Implicit residual smoothing is used to extend the stability range of the basic 
scheme. The residual is defmed by 



where m is the Runge-Kutta stage and D(m) is the total artificial dissipation at stage m. 
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Since the residual is only dependent on the neighboring cells in an explicit scheme, 
information at one cell may take many time steps to reach another cell. The speed at 
which this information would normally travel limits the maximum time step allowed by 
the scheme. Therefore, implicitly averaging the residuals makes the solution at any 
point partially dependent on the solution at all the other points. This increases the 
stability of the scheme and allows for larger time steps to be taken than would be 
allowed in the basic scheme. This method typically allows for a factor of two increase 
in the CFL number for the basic scheme, making the CFL number 7.5 for the scheme 
described. A further discussion of this method can be found in reference [60] . 

The main difficulty in finding a steady-state solution to the Navier-Stokes 
equations is the elimination of the unsteady acoustic waves generated during the 
initialization. The hardest waves to eliminate are the ones that have wavelengths that 
are the same order of magnitude as the grid spacing, or longer. To help eliminate these 
waves, multigrid methods have been developed. The purpose of multigrid is to change 
the number of grid points in an orderly manner so that waves with long wavelengths on 
the fine grids will have short wavelengths on the coarse grids. In this manner, waves 
that ordinarily would have taken many iterations to eliminate could then be eliminated 
in fewer iterations. 

Coarser meshes are obtained in the multigrid method by the elimination of every 
other mesh line in each coordinate direction. The solution is then transferred to the 
coarser mesh by a rule that conserves mass, momentum, and energy. Residuals are 
transferred to the coarser meshes by summation and are used as a forcing function to 
represent the fine grid solution on the coarse mesh. These strategies are discussed in 
more detail in references [61] and [62], This process is repeated until the coarsest mesh 
is reached. Then the corrections are transferred to the next finer mesh by bilinear 
interpolation. As stated in reference [59], the coarse grid corrections are smoothed 
using the implicit residual smoothing scheme discussed before, with constant 
coefficients, before they are passed to the next finer grid. Using this technique allows 
the multigrid method to be effective for a wider range of flows. 

A fixed W-cycle multigrid technique, with a full multigrid method, (FMG), is 
used to provide a well-conditioned starting solution for the finest mesh. These 
multigrid strategies are shown in figures 10 and 1 1 below. The implementation of the 
W-cycle in this code consists of a single sweep through the domain at each grid level. 
One complete multigrid cycle begins at point 1 in figure 10 and finishes at point 2. In 
the initial stages of the full multigrid method, shown in figure 1 1, 100 multigrid cycles 
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Figure 10-W-cycle Implementation for Multigrid 



Figure 1 1-FMG Implementation Method 


are completed at each of the coarse levels shown as points 1, 2, and 3. Point 4 signifies 
the beginning of the W-cycle for the calculation. By using this approach, acoustic 
waves with long wavelengths can be damped in the coarse levels to provide a good 
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initial solution to the finest grid level. 

Artificial viscosity is necessary in central differencing schemes for maintaining 
numerical stability. This viscosity will suppress high-frequency waves that are not 
damped by central differencing schemes and will dampen oscillations from capturing 
shock waves. Scalar dissipation models typically work well with most problems, but 
are not accurate enough for laminar flow problems. In this study, scalar dissipation is 
replaced with a matrix dissipation. This then gives the appropriate amount of artificial 
viscosity to each individual wave. This technique is fully described in reference [63] . 

The artificial dissipation technique used is similar to that originally developed 
by Jameson, Schmidt, and Turkel [57]. The total artificial viscosity consists of a blend 
of second and fourth differences in each coordinate direction. The purpose of the 
second difference term is to add an entropy-like condition to suppress oscillations near 
shock waves. Therefore, this term will be small in smooth regions of the flow field. 
The fourth difference is added to damp high-frequency waves and allow convergence to 
steady-state. The fourth difference, however, causes oscillations in the neighborhood of 
shocks. A shock detector, consisting of a pressure sensor, is used to turn off the fourth 
difference near shocks so only the second difference is operative. The pressure sensor 
used is given in reference [63], 

Finally, initial conditions and boundary conditions for the scheme must be 
addressed. The inflow and outflow boundary conditions employ point boundary 
conditions based on characteristic theory. For the subsonic inflow of the nozzles, four 
of the dependent variables are set from a one-dimensional analysis and the fifth is 
determined from the solution inside the domain. Zeroth-order extrapolation is done for 
the supersonic outflow boundary. The no-slip boundary condition is imposed at the 
solid walls. It is also assumed that the walls are adiabatic so the temperature is 
computed from the reduced energy equation 



(14) 


where t| is the direction normal to the surface. The wall pressure is calculated by 
setting the normal pressure gradient to zero at the wall. The initial conditions were set 
using simple one-dimensional area ratios. Using this technique, the pressure, density, 
temperature, and total velocity could be set at every point. The velocity components 
were then found by assuming that the velocity at every point was tangent to the local 
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grid. This technique gives a reasonable initial approximation to the final solution. 


2.2.3 Crossflow Determination 

After the Navier-Stokes calculations on the nozzles were completed, the 
boundary-layer crossflow needed to be calculated from the solution files. Since the 
flow is symmetric around both the centerplane and the diagonal, the crossflow only 
needed to be determined on one wall. The wall chosen for this is the one perpendicular 
to the z-direction. At each constant x and constant y location from the center of the 
nozzle to the edge of the boundary-layer on the wall that is perpendicular to the y- 
direction, z is varied from the wall to the edge of the boundary-layer to find the location 
of the maximum crossflow velocity (WmAX) and the distance between the wall and the 
location where the crossflow velocity becomes 10 percent of its maximum value (5 io). 

To calculate the crossflow, the unit vector in the crossflow direction must be 
determined. Since the crossflow must be perpendicular to the wall normal vector and 
the local inviscid flow, the dot product of the crossflow vector with both the wall unit 
normal and the local inviscid flow vector must be zero. This is given by 

qi n i + q2 n 2 + q3 n 3 = o (i 5 ) 

qi u e + fl2 v e + fl3 w e = 0 (16) 


where 

qm = component in the crossflow direction 
n m = component locally normal to the wall 
ue, v e , w e = component of the local edge velocity 

A discussion of the effects of the boundary-layer edge definition on the calculation of 
the crossflow can be found in the next chapter. To solve for qi, q2, and q3 uniquely, a 
normalization equation is needed. Assuming that most of the crossflow will be parallel 
to the wall and perpendicular to the x-direction, q2 is first chosen to be one. Solving the 
equations 15 and 16 with this assumption gives 

qi = ~ v <?~ q 3 w <? 

u e 


( 17 ) 
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q2 = i 


q 3 = 


v e ni-u e n2 

u e n 3“ w e n l 


(18) 

(19) 


The solutions to the above equations are then normalized to obtain a unit vector in the 
crossflow direction. 

This technique can also be used to determine the streamwise direction. Since 
this direction is perpendicular to the wall normal and the crossflow direction, the dot 
products are taken with these vectors and set to zero. Choosing s i to be one, since most 
of the streamwise profile is in the x-direction, the components of the streamwise 
direction becomes 


s i =i 

(20) 

S2 = -qi~q3 s 3 

(21) 

q2 


S3= amzS2EL 

(22) 


q2 n 3 " q3 n 2 

These equations are also normalized to obtain the unit vector in the streamwise 
direction. 

The wall normal needs to be found at every computational plane that makes up 
the z-wall. Since a vector is normal to a plane if it is normal to two lines in that plane 
that are not parallel, only two lines in each computational plane need to be known. The 
two lines that are chosen for this are 

ll=x^i + y^3 + z^k (23) 

1 2 =x T1 i + yi 1 j + z 11 k (24) 

Specifically, these two lines represent the vectors that pass through the points i,j,k and 
i+l,j,k and i,j,k and i,j+l,k respectively. The i counter changes in the streamwise 
direction and j in the direction parallel to the wall on which the crossflow is being 
calculated. Taking the dot product of both these lines with the wall unit normal gives 
two equations: 


x ^ n l+y^ n 2+z£n3=0 


(25) 
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xrjni + yii n 2 + z r| n 3 = 0 (26) 

Here again one more equation is needed to solve for ni, n2, and n 3 uniquely. 
Arbitrarily taking n3 equal to one, solving the set of equations, and normalizing them to 
obtain the wall unit normal gives 


n l = 


yrre-ypn 


(27) 


n 2 = 


n 3 = 


(y^ -y^) 2 +(*n*§ +(^y T1 -x T1 y^) 2 

~~ ( X T) Z ^ ~x^ z r)) 

(y^ -y^ z Ti) 2 +(x t 1 z ^-x ^ Zti ) 2 +(x^y n — x^)" 

~( x ^yri-^) 

(y n z ^ -y^ z n) 2 + (xn z ^ -x^n) 2 +(x^yti -Xqy^) 2 


(28) 


(29) 


where 

= x i+l,j,k ~ x i,j,k 

y^ = yi+i,j,k-yi,j,k 

= z i+l,j,k ~ z i,j,k 
X T) = x i,j+l,k “ x i,j,k 
yii = yi,j+i,k ~ yij.k 
z Tl = z i,j+l,k ~ z i,j,k 

x, y, z = stream wise and two normal-to-wall nozzle coordinates 
T|, C = streamwise and two normal-to-wall computational coordinates 

After the wall unit normal is known, the crossflow and streamwise directions 
can be calculated assuming the edge parameters are known. Finally, knowing the 
crossflow and streamwise directions, the magnitude of the crossflow and streamwise 
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profiles can be found at every point inside the boundary-layer by taking the dot 
products of the crossflow and streamwise directions with the local velocity vector. The 
magnitude of the crossflow can then be used to calculate the various crossflow 
Reynolds numbers used in this investigation. 

The procedure used for calculating the standard crossflow Reynolds number is 
straightforward. First, the crossflow magnitudes are searched to find the maximum 
crossflow velocity in any profile. After this point is found, the furthest distance from 
the wall, inside the boundary-layer, where the crossflow becomes 10 percent of the 
maximum is determined. Then, since the edge parameters are known, the crossflow 
Reynolds number can be calculated from equation 2. A region of s-shaped profiles 
existed in both nozzles that complicated the calculation of the parameter, however. If 
the above definition is used for all the parameters, the crossflow Reynolds number can 
still be calculated but a discontinuity occurs in the solution. This will be discussed 
further in the results for this parameter. Stability is determined at any point for this 
method by equation 30 developed by Godil and Bertelrud [40], 

R cftr =200|'l+^M*1 (30) 


where 

Rcffr Transition Crossflow Reynolds Number 
Me Local Edge Mach Number 


After the standard crossflow Reynolds number is known, Reed and Haynes 
crossflow Reynolds number can be calculated by the addition of the corrections for a 
cooled wall and compressibility as shown in equation 4. These corrections are 


H = 


Cfrio) 
C(&10)f T 

J j-\ 

0 vej 




(31) 
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L = 


K C ldJ 


3.279 + 1.721 


V 


1 w 

Tad 


|1 + A] + 0.664A 


(32) 


(5 + 2.385 A) 


where 


PrV 2 ( Y -l)M^ 

2 




Me Local edge Mach number 

FV Prandtl number 

C(8l0) Distance along wall normal from wall to location of 5 io 

T w Local wall temperature 

T/T e Local temperature divided by local edge temperature 

H Correction for compressibility effects 

L Correction for wall cooling effects 
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Since the wall and edge parameters are known at any point, the correction for wall 
cooling is easily calculated. The integral of T/T e is calculated using a trapezoidal 
integration method across the boundary-layer, along the wall normal. The 
compressibility correction can then be calculated by dividing the distance from the wall 
to the edge of the boundary-layer by the integrated boundary-layer temperature. 
Multiplying these two corrections with the standard crossflow Reynolds number gives 
the Reed and Haynes crossflow Reynolds number. Stability for this method is 
determined by the equation 


R = 


R cfR&H 

Wmax> 

/U e 


(34) 


where W MAX/Ue is in percentages. Stability is then predicted for R less than or equal 
to 44.0 for quiet flows. This method is limited, however, to the range 2% < WMAXAJe 
<8% corresponding to the range of data that was investigated when defining the 
parameter R. 

The crossflow Reynolds number developed by Beasley can be calculated by 
integrating the crossflow magnitude from the wall to the edge of the boundary- layer, 
along the wall normal, using a trapezoidal integration technique. This result is then 
divided by the kinematic viscosity at the edge as in equation 3 to get the solution. Also 
of interest is the streamwise shape factor (H s hp) for stability analysis. First, both the 
displacement thickness (6*) and the momentum thickness (0) are found by integrating 
the streamwise velocity profile across the boundary-layer using the compressible 
formula [21]. The integration is again done using a trapezoidal integration technique. 
Finally, the shape factor is calculated by dividing the displacement thickness by the 
momentum thickness. The stability limits of the solution is then set by equation 35 
which is valid for the range 2.3 < H s hp < 2.7. 


Cl = 


300 f -i 
tan 


71 


0.106 


Khp -2.3) 


2.05 


(35) 


A computer program has been written following these procedures to determine 
the crossflow in both Mach 2.4 nozzles. A listing of the program can be found in 
appendix B, a brief description of the methods follows. The program reads in two 
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binary files consisting of the geometry file created by the grid generation program and 
the solution file of dependent variables created by the Navier-Stokes solver. The 
program writes four files. First, DIFFP is the difference in pressure between the comer 
and the center of a wall. This data is used to compare with the original design data. 
Second is the file titled CROSSFLOW and contains the data for contour plots of the 
crossflow Reynolds numbers and the stability criterion. Third is a file entitled 
CORNER which creates a comerplane cut of the data for plotting Mach contours. 
Finally, CROSSFLOWMAG contains line plots of the crossflow magnitude in the 
boundary-layer at various locations. 

After the data is read in and manipulated, the edge of the boundary- layer on the 
wall perpendicular to the y-direction is found on the centerline. This is the outer limit 
for the calculation of the crossflow so the comer would not be included in the 
calculation. Then, at each point, the wall normal is calculated and a check is made to 
make sure it is inside the computed domain. The program then marches along the wall 
normal and saves the velocity every time the normal vector crosses a grid line. This 
procedure continues until the boundary-layer edge is crossed. After the edge is found, 
the stream wise and crossflow directions can be calculated. The program again marches 
through the boundary-layer to calculate the streamwise and crossflow magnitudes, find 
WmAX and 5 10* and integrate to find the compressibility correction, displacement 
thickness, and momentum thickness. Finally, the crossflow Reynolds numbers and the 
stability criterion are calculated and the program ends. The data is then plotted and 
analyzed to get an estimate of the boundary-layer stability. 
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CHAPTER m RESULTS AND DISCUSSION 


The flow fields in both Mach 2.4 nozzles were simulated with the 3D Navier- 
Stokes solver LARCK. The approximate memory and run-time requirements to 
converge the solutions by four orders of magnitude, starting from the first FMG level, 
on NASA Langley's CRAY-YMP are given below: 


Table #2 - Computational Requirements 


Grid 

Hours 

Memory 

(Megawords) 

137x41x41 

12 

15 

201x65x65 

30 

48 

305x97x97 

142 

158 


Increasing the number of grid points has a linear effect on the amount of memory 
required. The time necessary for convergence, however, is not linear. This is most 
likely due to the increase in cell aspect ratios in the finer grids. The convergence 
behavior for these solutions is shown in figures 12-15. Figure 12 is the convergence 
history for the short Mach 2.4 nozzle using the coarse grid, figure 13 is the medium 
grid, and figure 14 is the fine grid. Finally, figure 15 is the convergence history for the 
long nozzle using the 201x65x65 grid. The discontinuities in each of these figures 
occurs with the addition of another level to the FMG cycle. When a new level is added, 
an interpolation of the solution to the new level must be performed. This interpolation 
introduces round-off error into the solution, causing the jump in the residual. The 



















LOG 10 (Error) „ LOG 10 (Error) 
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round-off error is, however, eliminated in approximately 50 multigrid cycles after the 
introduction of a new level. Also note that four orders of magnitude reduction in the 
residual was not achieved in the fine grid calculation. The poor convergence rate for 
this grid made it too expensive to complete. This was most likely due to the high aspect 
ratios of the grid, along with the possibility of poor damping at the inlet of the nozzle. 
It was felt that the 3 to 3.5 orders of magnitude that was achieved made the solution 
adequate for comparison with the other grids. 


11 Boundarv-Laver Ed ye Definition 

The first step in analyzing the data was to determine how the boundary-layer 
edge definition effected the calculation of the boundary-layer parameters. Since this 
definition effects the calculation of the crossflow Reynolds numbers, the definition is 
important. The boundary-layer edge definition analysis was done using the results from 
the medium resolution grid on the short Mach 2.4 nozzle. 

Since the streamwise velocity in the 'in viscid' flow region need not be constant 
at any streamwise location, the edge location is relatively arbitrary in a Navier-Stokes 
calculation. However, in shock-free nozzle flow, the total pressure must remain 
constant outside the boundary-layer. Therefore, this parameter can be used to determine 
the boundary-layer edge. The question then becomes how much of a total pressure drop 
must exist before the edge is reached. For a three-dimensional nozzle with no 
experimental results to verify analysis (like those presented here), this question becomes 
even more difficult. 

To investigate how the velocity and total pressure change through the boundary - 
layer, the u-component of the velocity was plotted along with the total pressure divided 
by the frees tream total pressure near the wall. This can be seen in figure 16. Note that 
these values are taken from the centerplane of the nozzle at the exit. At this location, 
the wall normal is perpendicular to the x-direction and therefore, the u-component of 
the velocity should be equal to the total velocity in the streamwise direction. Therefore, 
by using this location, the analysis is greatly simplified. In this figure, the u-velocity 
appears to reach the edge at a y-location between 0.107m and 0.109m. These locations 
correspond to the values of Po/P 0 FS of between 1.0 and 0.9, respectively. 

Two edge locations, P 0 = 0.90 PoFS and 0.98 PoFS. were then chosen to further 
analyze the crossflow Reynolds number data. The data for WMAX/Ue in percentages 
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Y(m) 

Figure 16-Total Pressure Ratio and U- Velocity Component in the 
Boundary-Layer Region 


along half of the z-wall for both edge definitions is plotted in figure 17. Note the line 
approximately 0.4m downstream of the throat that is unlabeled. This line is the 
interpolated location where the crossflow changes direction from positive, or towards 
the comer, to negative, or away from the comer. The direction change causes a jump in 
both the standard and the Reed and Haynes crossflow Reynolds numbers, as can be seen 
in figures 18 and 19 respectively. Further analysis of the flow characteristics in this 
region are presented in the next section. 

Investigation of figure 17 shows that the edge definition in the nozzle throat 
region has a significant affect on the calculation of the maximum crossflow, while little 
difference is seen near the exit. The differences in the throat are most likely due to the 
differences in the edge velocity for each definition. Since the boundary-layer at the 
throat is relatively thin, the stream wise velocity changes rapidly over a small distance. 
Therefore, large differences can occur in the calculated edge velocity for each edge 
definition. In the exit, however, the boundary-layer is much thicker, and therefore, 
should have much smaller differences between edge velocities for each edge definition. 



Figure 17-WMAX/Ue in Percentages for T 


0.90 P 
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Figure 18 -Standard Crossflow Reynolds Number for Two Boundary-Layer Edge Definition 
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Figure 19-Reed and Haynes Crossflow Reynolds Number for Two Boundary-Layer Edge Definitions 
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Figure 20-Crossflow Reynolds Number by Beasley for Two Boundary-Layer Edge Definitions 
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Figure 21-Streamwise Shape Factor for Two Boundary-Layer Edge Definitions 
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This leads to smaller differences in the calculated crossflow, which was observed in 
figure 17. 

The various crossflow Reynolds numbers have been plotted in figures 18-20 for 
both of the boundary-layer edge definitions. Figure 18 shows the standard crossflow 
Reynolds number. As in figure 17, the standard crossflow Reynolds number is 
significantly affected by the change of edge definitions the in the throat region. 
However, near the exit, change in the crossflow Reynolds number between the two 
definitions is much smaller than in the throat. Assuming that the edge definition does 
not affect the kinematic viscosity very much, the changes seen in the throat are due to 
the changes in the maximum crossflow velocity. The differences in W MAX affect 8io 
in this calculation as well. Figure 19 shows Reed and Haynes crossflow Reynolds 
number for both edge definitions. This figure shows the same trends that were seen in 
the standard crossflow Reynolds number. Therefore, the edge definition does not seem 
to significantly affect the corrections for compressibility or wall cooling in Reed and 
Haynes correlation. Finally, figure 20 shows the crossflow Reynolds number developed 
by Beasley, shown in equation 3, for both edge definitions. Since this parameter 
depends on the integral of the crossflow velocity profile, the effect of the edge 
definition on the crossflow direction can be better estimated from these results than 
from the previous two crossflow Reynolds numbers. This is assuming, however, that 
the edge definition does not affect the kinematic viscosity very much at the boundary - 
layer edge. Near the throat, Beasley's crossflow Reynolds number seems to be severely 
affected when the edge definition is changed. However, near the exit, the change in the 
edge definition affects the results much less than in the throat region. There appears to 
be more of a change in this parameter near the exit than in the other crossflow Reynolds 
numbers, however. This signifies that the change in the edge definition may affect the 
crossflow calculation more than was apparent from the previous results. 

The streamwise shape factor is plotted for both edge definitions in figure 2 1 . 
This figure shows that only in the throat region does the edge definition seem to affect 
the streamwise shape factor. Downstream of where the streamwise shape factor is 3.0, 
the two calculations are almost identical. However, the throat region is the region of 
interest for this parameter since the Cl prediction criteria is only valid for 2.3 < H s hp < 
2.7. 

Obviously, the boundary-layer edge definition chosen for a Navier-Stokes 
calculation can significantly change crossflow Reynolds number results. Although it 
has been shown that the boundary-layer edge definition is important in calculating the 


57 


crossflow Reynolds number and estimating transition, the results presented do not help 
determine which edge definition is the most suitable. King [17] has previously 
performed a three-dimensional Navier-Stokes calculation for a two-dimensional nozzle. 
This investigation found that using equation 36, below, gave excellent results compared 
with results from a two-dimensional boundary-layer calculation. Results were also in 
agreement with experiments for 6* but were not as accurate for 0. 

Po edge = °- 98po max +0 02p <>wall (36) 


where 

Poedge 

Pomax 

Powall 


Total pressure at the boundary-layer edge 
Maximum calculated total pressure 
Total pressure at the wall 


Referring to figure 16, it can be seen that the total pressure at the wall is only about 8% 
of the freestream value. Substituting this into the above equation and approximating the 
maximum total pressure as the freestream total pressure, the edge location is found to be 
approximately where the total pressure is 98% of the freestream value. Therefore, this 
is the definition that will be used for this investigation. 


3.2 Short Mach 2.4 Nozzle Results 

The results for the short Mach 2.4 nozzle are separated into two sections. First, 
a grid resolution study was performed to ensure that the calculated results are 
independent of the grid spacing. The second section then presents crossflow and 
transition results as well as other flow properties of interest. 


3.2.1 Grid Resolution 

The Mach number distributions along the centerline for the short Mach 2.4 
nozzle for all grids and for the initial axisymmetric design are shown in figure 22. 
Good agreement exists between all solutions and the initial design. Therefore, in the 
streamwise direction, the coarse grid distribution should be all that is required to resolve 




Mach Number 

Figure 23-Mach Number Distribution at the Exit, in the Centerplane, All Grids 




Li! 


0.1005 
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201x65x65 Grid 
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Mach Number 

Figure 24- Mach Number Distribution Near the Wall in the Centerplane, 

at the Exit, All Grids 
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Figure 25-Pressure Difference Between the Corner and Centerplane, 
All Grids and In viscid Design 
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this portion of the flow. Figure 23 and 24 show the Mach number distribution at the 
centerplane, at the exit. In figure 23, it appears that there is good agreement between all 
solutions. However, figure 24 shows that in the boundary- layer region there is a 
discrepancy in the solutions. It can only be assumed that the fine grid gives the most 
accurate results, since further verification is not currently possible. Finally, in figure 
25, the difference in pressure between the comer and the centerplane is shown for all 
three grids and the initial inviscid design. It can be seen that there is a good agreement 
between all three grids. However, the Navier-Stokes solutions do not agree very well 
with the inviscid design. This is most likely due to boundary-layer effects, such as 
crossflow and comer flow. 

The displacement and momentum thicknesses were also computed along the 
centerplane, which are shown in figures 26 and 27. Both of these plots show a general 
agreement between the Navier-Stokes solutions up to approximately 0.5 meters 
downstream of the throat. At this point they diverge to a maximum difference at the 
exit. The table below gives the percentage difference from the finest grid for all the 
solutions and the two-dimensional boundary-layer results, at the exit: 


Table #3 - Percent Difference in 8* and 0 from Finest Grid 



5* 

0 


28.30 


137x41x41 

26.42 

15.47 

it) 1x65x65 

9.43 

7.14 

mmsRsm 

0.0 



Again, it can only be assumed that the fine grid solution is the most accurate 
calculation. One possible reason for the differences could be variations in the resolution 
with which the crossflow is captured. To further investigate this, the streamwise and 
crossflow profiles were plotted at the exit, half-way between the centerplane and comer. 
These are shown in figures 28 and 29. Figure 28 shows the streamwise profiles. It can 
be seen that increasing the resolution makes the profile more shallow. The decrease in 
velocity in the streamwise direction is accompanied by an increase in crossflow as is 
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Figure 26-Displacement Thickness in the Centerplane, All Grids and 
Two-Dimensional Calculation 



X(m) 

Figure 27-Momentum Thickness in the Centerplane, All Grids and 
Two-Dimensional Calculation 
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Figure 28-Streamwise Profiles at the Exit, Approximately Halfway 
Between the Centerplane and the Comer, All Grids 



Figure 29-Crossflow Profiles at the Exit, Approximately Halfway 
Between the Centerplane and the Comer, All Grids 
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shown in figure 29. This increase in crossflow will also increase the mass and 
momentum transfer from this location towards the centerline. Therefore, as the 
resolution is increased, the increase in mass and momentum transfer towards the 
centerline causes the higher displacement and momentum thicknesses seen in figures 26 
and 27. It should also be noted that the differences in figures 28 and 29 correspond to 
the differences in Mach number near the wall that was shown in figure 24. 

Two other points should be noted for figures 26 and 27. First, there is a 
relatively good agreement between all the grids and the two-dimensional results up to 
the throat region for both 5* and 0 . Using figure 30, it can be seen that crossflow plays 
a small role up to the throat region, and therefore, transfer of mass and momentum into 
the centerplane boundary-layer should be small. Therefore, the two-dimensional results 
should be in line with the three-dimensional results up to the throat. This signifies that 
using P 0 = 0.98 PoFS was a good choice for the boundary-layer edge. Second, after the 
throat, the slopes of the lines for the computed 6* and 0 are smaller than the two- 
dimensional results up to x=0.5m. Since the crossflow is toward the comer in this 
region, mass and momentum are transferred away from the centerplane boundary-layer. 
This results in smaller growth rates for 6* and 0 than the two-dimensional results 
predict. At approximately x=0.5m downstream of the throat the slopes of these lines 
becomes greater than the two-dimensional results. Again using figure 30, it can be seen 
that this is the approximate location where the change in crossflow direction from 
towards the comer to away from the comer begins to effect the centerline boundary - 
layer. Opposite from before, mass and momentum are now being transferred into the 
centerline boundary-layer, causing larger growth rates for 8* and 0 than predicted by 
the two-dimensional results. 

Figure 30 shows the maximum crossflow velocity (WmAX) divided by the local 
edge velocity (Ue) in percentages. In each plot there is an unmarked line approximately 
half way between the throat and the exit. This is the interpolated location where the 
crossflow velocity changes sign from positive (towards the comer) to negative (away 
from the comer). In this region s-shaped profiles exist, which create a discontinuity in 
the crossflow Reynolds numbers. This will be discussed in more detail in the next 
section. It can also be seen that increasing the resolution increases the maximum 
crossflow velocity, assuming that the edge velocity is approximately constant at any 
given location. The medium and fine grids have approximately the same results with 
the major difference being the larger area in the fine grid calculation where W MAX Ale 
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Figure 30- Wj 
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Figure 31 -WmaxAJ e in Percentages Near the Corner, at the Throat 


is -6%. The coarse grid, however, only reaches a maximum crossflow of - 5 % near the 
exit. 

One further note should be made with respect to the throat data in figure 30. As 
can be seen in this figure, the crossflow appears to increase in magnitude into the 
comer. This should not occur, however, since the wall perpendicular to the crossflow 
should cause the crossflow to go to zero in the comer. To further investigate this, the 
contours have been magnified near the comer, in the throat region, for the medium 
resolution grid. This is shown in figure 31. This figure shows that the crossflow 
velocity does in fact decrease to zero and even reverses slightly as the corner is 
approached. 

The various crossflow Reynolds numbers and transition estimation methods are 
plotted in figures 32-36. The meaning of these will be discussed in the next section and 
are presented here only to determine grid effectiveness in their calculation. In figure 32, 
it appears that increasing the resolution only increases the distinctiveness of the jump 
location created by the discontinuity in 5 10- Near the exit, the contour corresponding to 
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a crossflow Reynolds number of 2000 becomes more clearly defined with increasing 
resolution. This is most likely due to the increased resolution of WMAX seen in figure 
30. The same trends are also seen for Reed and Haynes crossflow Reynolds number in 
figure 33. However, the values in the throat region change more rapidly than the 
standard crossflow Reynolds number does, and the contours also become more defined 
with increasing resolution. The changes near the throat for RcfR&H must be due to 
changes in the compressibility factor since similar trends are not seen in figure 32. 
Finally, the crossflow Reynolds number developed by Beasley is seen in figure 34. 
This figure also shows the same trends that were seen in the standard crossflow 
Reynolds number. However, since this variable depends on the integral of the 
crossflow profile, the profiles seem less dependent on the resolution. Note that when 
the exit region is reached, the profiles becomes fuller with increasing resolution. These 
fuller profiles were seen in figure 29. 

The transition estimation parameter R is shown in figure 35 for all grid 
resolutions. It can be seen in this figure that the medium and fine grids generally agree 
with the increase in resolution, only making the jump location more defined. However, 
the coarse grid does a poor job in the calculation of this parameter with respect to the 
other two resolutions. The streamwise shape factor is shown in figure 36 for all grid 
resolutions. There are discernible differences between all grids and especially between 
the coarse grid solution and the medium and fine grid solutions. These differences in 
streamwise shape factor are not surprising considering the differences in displacement 
thickness and momentum thickness shown in figures 26 and 27. 

Since the coarse grid cannot accurately represent all parameters necessary for 
transition estimations, it cannot be used for these studies. The fine grid provides what 
are assumed to be the most accurate results, but is too costly in both memory 
requirements and run-time to allow repeated use. The medium grid will be used for the 
studies presented in this paper since the medium grid results are similar to the fine grid 
results. Clearly, the resolution is insufficient for stability calculations. However, the 
resolution seems more than sufficient for preliminary investigation using the crossflow 
Reynolds numbers, since they are fairly consistent. 
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Figure 32-Standard Crossflow Reynolds Number, All Grids 
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Figure 33-Recd and Haynes Crossflow Reynolds Number, All Grids 
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Figure 34-Crossflow Reynolds Number of Beasley, All Grids 
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Figure 35-Reed and Haynes Transition Estimator, R, All Grids 
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Figure 36-Streamwise Shape Factor, All Grids 
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3.2.2 Results 

After the Navier-Stokes solution was complete, the Mach contours were plotted 
in both the centerplane and the comerplane. This can be seen in figure 37. It can be 
seen that the test rhombus is fairly uniform with a maximum deviation of +0.017 from 
the test Mach number of 2.4. The Mach number in the exit plane is then plotted in 
figure 38. This plot shows that the Mach number at this location is also very uniform 
with the maximum deviations being ±0.0025 from the design Mach number of 2.4. A 
small comer flow region should also be noted which contains about 3% of the total wall 
length. Figure 39 shows the v and w velocity vectors in this comer region. As can be 
seen, two counter-rotating vortices are present. This is most likely the early stages of 
the large counter-rotating vortices that were present in the hypersonic rapid-expansion 
nozzles investigated in reference 47. It is possible that these vortices will cause 
transition near the comer of the nozzle. However, because of their relatively small size 
and the fact that the maximum velocity is only around 1% of the u-component of the 
velocity, these vortices are not likely to dominate the transition. 

Flow in the comer of these nozzles is complex and three-dimensional, and the 
crossflow Reynolds number techniques are not applicable. Therefore, the physical 
limits for using the crossflow Reynolds numbers must be set before a quantitative 
analysis can be completed. To accomplish this, the boundary-layer thickness has been 
plotted from the centerplane into the comer at both the throat and the exit. These are 
shown in figures 40 and 42. Figure 40 shows that in the nozzle throat, the boundary - 
layer remains fairly constant up to about 0.0013m from the wall. The thickness 
suddenly jumps after this point. A close-up of just the jump region is shown in figure 
41. This plot shows that there is a distinct increase in the boundary-layer thickness at 
y=0.0014m from the comer. At this point, the comer flow region is entered and the 
crossflow data is no longer useful. To further see this, the boundary-layer in the exit 
plane is shown in figure 42. The figure shows that the boundary-layer at the exit grows 
until a jump again occurs where the boundary-layer becomes much thicker. Figure 43 
shows a close-up of this region where it can be seen that the boundary-layer thickness 
increases by approximately an order of magnitude when the comer region is entered. 
To show that this jump is created by the flow in the comer, figure 44 shows total 
pressure contours near the comer region. As can be seen, the boundary-layer is found 
to be approximately constant up to the point labeled "jump location". After this point, 
the comer flow properties take-over and the crossflow Reynolds numbers are no longer 
valid. Note that the point labeled "jump location" is at approximately 0.0035m from the 
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Figure 37-Mach Contours in the Centerplane and Comerplane 
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Figure 38-Mach Contours at the Exit 
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Figure 39-V and W Velocity Vectors in the Comer at the Exit 
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Figure 44-Total Pressure Contours at the Exit, in the Comer 


comer, the same location where the discontinuity occurs in the data of figure 43. The 
limits for the crossflow calculations in the short Mach 2.4 nozzle were taken to be 
0.0014m from the comer in the throat and 0.0035m from the comer at the exit. 

As mentioned in the previous section, there exists a region in these nozzles 
where the crossflow changes sign from positive to negative. Figures 45 and 46 show 
streamwise and crossflow profiles from x=-0.927 into the nozzle exit, along with the 
line along which the profiles were taken. It can be seen in figure 45 that the streamwise 
profiles generally increase in magnitude until approximately 0.4 meters downstream of 
the throat. At this location, the streamwise velocity becomes approximately constant to 
the exit. In figure 46, it can be seen that the crossflow is positive (towards the comer) 
until somewhere between x=0.3 and x=0.4 meters downstream of the throat, where the 
profiles undergo a change in direction. To better see this, figure 47 shows the profiles 
through this region. Note that the profiles are plotted at intervals of Ax = 0.00658m. 
Since both Ref and RcfR&H are dependent on the maximum value of the crossflow and 
the location where the crossflow decreases to 10 percent of the maximum value (5io), 
the calculation of these parameters in this region is difficult. The following definition 
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has been chosen for this region: after WmAX is found, the location for 5io is taken to 
be the furthest point from the wall where the magnitude of the crossflow is 10 percent 
of the maximum. As can be seen in figure 47, the initial crossflow profile shown is 
entirely in one direction. Then, moving downstream, the crossflow begins to change 
direction near the wall while continuing to slow down near the edge. As this happens, 
WmAX decreases to a minimum and then begins to increase again. Since 5io is 
defined to be the furthest distance away from the wall, this parameter will continue to 
increase until the maximum value is reached at point A in figure 47. After this point, 
the crossflow in the top half of the s-profile is always less than 10 percent of WmAX 
Therefore, the distance 5 io will suddenly jump to point B shown on figure 47. This 
then causes the discontinuity seen in the crossflow Reynolds numbers. 

The crossflow Reynolds number contours as well as the transition estimation 
criteria are shown in figures 32-36 along half of any nozzle wall. The standard 
crossflow Reynolds number, plotted in figure 32, shows two regions where transition 
may occur. The first, and most significant, is near the exit, where R c f increases from 
approximately 150 to more than 2000 over the range of about 0.16m. The maximum 
crossflow Reynolds number attained is almost five times as large as the estimated 
critical value of 430. It is thus likely that transition will occur before the exit. The 
second region is that near the throat. Although a maximum of slightly more than 500 is 
attained in this region, the estimated transitional value of 240 is approximately 1/2 of 
the calculated value. Therefore, it is likely that transition will occur near the nozzle 
throat, according to this method. 

Figures 33 and 35 show the Reed and Haynes crossflow Reynolds number and 
the transition estimator R. It can be seen in figure 33 that, like R c f, there is a region 
near the exit where the values change from about 250 to about 1250 over a region of 
approximately 0.16m. A region of lesser values also exists in the throat, like those seen 
previously. The values for R shown in figure 35 give the same general conclusions as 
the standard crossflow Reynolds number. After the observed jump region between 
x=0.4m and x=0.5m, the maximum values for R increase to more than 250. This is 
more than 5.5 times greater than the value of 44 estimated for transition on spinning and 
yawed cones [39]. Again, transition appears likely before the exit. Near the throat, the 
values for R become 2 to 3 times as great as the estimated transitional value. However, 
referring back to figure 30, this method can only be used in the small region where 
WMAX/Ue is greater than 2% to be within the lower limits for this correlation. Still, 
this method predicts that there is a good possibility that the flow will transition in the 
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Figure 47-Crossflow Profiles Through S -Shaped Region 


throat region. 

Finally, figures 34 and 36 show R§2 and the streamwise shape factor 
respectively. Figure 34 shows the same trends that were evident in the other two 
crossflow Reynolds numbers. A comparison of figure 34 with figure 33 reveals that the 
crossflow Reynolds number developed by Beasley produces approximately the same 
results as Reed and Haynes crossflow Reynolds number. Examination of the 
streamwise shape factor shows that this method estimates transition as soon as R 52 

reaches a value of 150. It is evident that transition will occur in the throat for this 
method. Since the limit on the streamwise shape factor is 2.3 < H s hp < 2.7, this method 
cannot be applied anywhere downstream of the throat. Therefore, no conclusions can 
be drawn from this method downstream of the throat area. 

Even taking into consideration the large discrepancies that may exist in these 
methods, there is little doubt that crossflow induced transition will occur in this nozzle. 
The table below shows the estimated transition location for each method both near the 
throat and near the exit. The transition estimates are those presented in the preceding 
paragraphs where the locations used are the first streamwise locations where any point 
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is greater than the transition estimate. 


Table 4-Transition Estimates for the Short Mach 2.4 Nozzle 


Method 

Location Near 
Throat 

Location Near 
Exit 

Ref 

X=-0.045m 

X=0.42m 

RcfR&H 

X=-0.1m 

X=0.4m 

R62 

X=-0.05m 

Inapplicable 


Clearly, the first and third method approximate transition at the same location near the 
throat while Reed and Haynes is not as close. However, Reed and Haynes estimation 
was taken from the lower end of the range on WMAXAJe and therefore, may not be 
reliable. Near the exit, the first and second method estimate transition in about the same 
location while nothing can be said about the third technique. It should be noted that the 
transition locations reported first occur near the comer for all cases. 

Even though transition is likely in the throat region, the likelihood is greater that 
transition will occur near the nozzle exit. If only the results near the exit are considered, 
a simple preliminary analysis of quiet-test length for the nozzle can be performed. 
Figure 37 shows the Mach contours for the short nozzle in both the centerplane and 
comerplane. From this it can be seen that the test rhombus starts at approximately 
x=0.52m. If a characteristic is drawn beginning from the reported transition location 
above, in the comerplane plot in figure 37, whose slope is equal to the average Mach 
number though which the characteristic travels, the quiet-test length was found to be 
approximately 0.17m. This then gives a quiet-test Reynolds number (RAx) of 12 x 

107. This is comparable to other quiet- flow nozzles in use today [23]. 
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Figure 48-Mach Number Along the Centerline, Navier-Stokes Solution and 

Axisymmetric Design 


3.3 Long Nozzle Results and Comparisons 
A simulation was then performed on the long nozzle case using a 201x65x65 
grid with the same stretching used in the previous case. Since the exit dimensions and 
throat dimensions are the same, using the same number of points and the same 
stretching should give reasonably accurate results normal to the streamwise direction. 
However, since this nozzle is almost twice as long as the previous one, using 201 points 
in the streamwise direction yields almost twice as much distance between grid points. 
This spacing is only 17 percent larger than the coarse grid in the previous example, 
however, which had good accuracy in the streamwise direction. Therefore, it will be 
assumed that this distribution is accurate enough for this investigation. To verify this 
assumption, the centerline Mach number is plotted for both the Navier-Stokes solution 
and the original design in figure 48. It can be seen that there is good agreement 
between the original design and the final solution signifying that this spacing is 
adequate. A further grid refinement study will be needed if further analysis is to be 
done on this design. 
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Figure 49- Mach Contours in the Centerplane and Comerplane 
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Figure 52-Pressure Difference Between the Centerplane and Corner, Navier-Stokes 

Results and Inviscid Design 


After completing the Navier-Stokes calculation on this nozzle, the Mach 
contours were again plotted in both the centerplane and cornerplane. This can be seen 
in figure 49. This figure shows that the maximum deviation in the test rhombus is about 
+0.008 from the design Mach number of 2.4. This is about 50% less than in the short 
design signifying a more even flow in the test section. The Mach contours were then 
plotted at the nozzle exit plane, shown in figure 50, to investigate the uniformity. The 
maximum deviation at this location is +0.002 in Mach number which is in the same 
range as was seen in the short nozzle. Finally, figure 51 shows the v and w velocity 
vectors at the nozzle exit near the comer. Again, counter-rotating vortices are present 
like those seen in the short nozzle. However, the reference vector shows that these 
vortices are smaller in magnitude than those in the short nozzle. It is reasonable to 
conclude that a longer nozzle design will have smaller vortices near the comer. Also 
note that these vortices cover a region of 3% of the length of the wall and are therefore 
are not expected to dominate the boundary-layer transition. 

The difference in pressure between the comer and centerplane on the wall is 
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shown in figure 52 for both the inviscid design and the Navier-Stokes solution. As in 
figure 25, a large dip is seen in the nozzle throat, after which, the difference goes to 
about zero. Also note that the Navier-Stokes results seem to agree better with the 
inviscid design than in the short nozzle. There are significant differences between the 
two however, especially after the throat, where the inviscid design predicts another 
significant change not seen in the Navier-Stokes results. The reason for this dip is 
unknown, but could be caused by the Mach number differences at the beginning of the 
radial flow region. Further comparisons will be presented in the next section where an 
attempt is made to analyze the reasons for the large gradient in the throat. 

Figure 53 shows the values for WMAXAJe in percentages for the long nozzle. 
Comparison with the short nozzle results show that there are no significant differences 
in the throat region. However, near the exit there is a reduction in the maximum 
crossflow velocity of about 2 percent of the edge velocity, assuming that the edge 
velocities in the two nozzles are approximately the same at relative streamwise 
locations. Also evident is the more gradual build-up of crossflow in the exit region, as 
signified by the larger spacing between contours. Further investigation of this is shown 
in figures 54 and 55 where the streamwise and crossflow profiles are plotted, along with 
the line from which they are taken. Comparing figure 54 and figure 45, it can be seen 
that the profiles are approximately the same in both nozzles up to the throat. After the 
throat, there is a more gradual increase in the streamwise velocity in the long nozzle, as 
is expected in a slower expansion nozzle. Comparing figure 55 with figure 46, it is 
evident that the contours in the throat region are similar while those near the exit are 
approximately 2/3 the magnitude of those in the short nozzle. This can be expected, 
considering the data seen in the WMAX/Ue contours in figure 53 and 30. Finally, 
figure 56 shows the progression of the crossflow through the s-shaped region. Note that 
the profiles are taken at intervals of Ax = 0.02941m up to the dashed line and Ax = 
0.0147m after the dashed line. Here again, compared to the short nozzle, these contours 
are much more gradual in their change from one direction to the other. 

Before analyzing the crossflow Reynolds number data, the limits for which the 
data is useful must be known. The boundary-layer thickness is shown for the long 
nozzle at the throat in figure 40. It can be seen here that the boundary- layer thickness in 
the long nozzle is slightly thicker than in the short nozzle. This is expected since the 
long nozzle has a slower convergence to the throat, giving the boundary-layer more 
time to develop. This figure the same discontinuity in the long nozzle that was seen in 
the short nozzle. Figure 41 shows a closer view of the discontinuity region. From this 
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Figure 54-Streamwise Profiles and Line Along Which They are Taken 
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Figure 55-Crossflow Profiles and Line Along Which They are Taken 
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Figure 56-Crossflow Profiles Through S-Shaped Region 


figure, the comer flow begins at y=0.0014m from the comer. This is the same location 
as in the short nozzle. The boundary-layer thickness at the exit is shown in figure 42. It 
can be seen that the long nozzle boundary-layer is again thicker than in the short nozzle. 
However, this is expected since the long nozzle boundary-layer has almost twice the 
length in which to develop. The discontinuity seen in figure 42 is expanded in figure 
43. This figure shows that the comer region in the long nozzle begins slight closer to 
the comer itself, at y=0. 0020m from the comer opposed to 0.0035m seen in the short 
nozzle. 

Figures 57-59 show R c f, RcfR&H. and R 52 respectively. Comparing these to 
the short nozzle results, in figures 32-34 respectively, shows that the results are 
comparable in the throat region while the maximums are all about 25 percent smaller 
near the exit. This is consistent with the reduction seen in WMAXAJe for the long 
nozzle when compared with the short nozzle results. A similarity is again noticed 
between Reed and Haynes crossflow Reynolds number and Beasleys crossflow 
Reynolds number in figures 58 and 59 respectively. 
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Figure 60-Reed and Haynes Transition Estimator, R 
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Since the long and short nozzles are designed for the same Mach number, the 
transition crossflow Reynolds numbers are again 240 near the throat and 430 at the exit. 
Figure 57 shows that there is a region in the throat where the computed crossflow 
Reynolds number is greater than 500, or more than two times the estimated transition 
criterion. In the exit, the computed crossflow Reynolds number is greater than 1500, or 
larger than 3.5 times the estimated transition. The transition parameter R given by Reed 
and Haynes is shown in figure 60. An interesting point here is that the values for this 
parameter are larger than in the short nozzle while the crossflow Reynolds numbers are 
smaller. Although most of the data for the R parameter is greater than the predicted 
transition value of 44.0, most of the data is out of the suggested range of 2% < 
WMAX/Ue < 8%. Taking the data from the prescribed region, the transition parameter 
in the throat is 100, more than two times the estimated transition of 44.0. Near the exit, 
the transition parameter becomes a maximum of 300, which is almost 7 times the 
estimated transition value. Finally, the streamwise shape factor is shown in figure 61. 
Using the Cl criterion, transition should again occur as soon as R 52 becomes larger 

than 150. This occurs right in the nozzle throat and near the exit. However, the 
streamwise shape factor is only in the range of 2.3 < Hshp < 2.7 in the throat region. 
Therefore, transition can only be estimated in the nozzle throat, after which, no further 
conclusions can be drawn. The first occurrences of transition estimated by these 
various methods is shown in table 5 below 


Table 5-Transition Estimates for the Long Mach 2.4 Nozzle 


Method 

Location Near 
Throat 

Location Near 
Exit 

Ref 

X=-0.045m 

X=0.975m 

RcfR&H 

X=-0.1m 

X=1.0m 

R62 

X=-0.075m 



None of the methods have comparable results for transition near the throat. It is 
unclear why the differences occur in these methods near the throat. Near the exit, the 
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first and second method again estimate transition in about the same location while 
nothing can be said about the third technique. However, the two methods are not as 
close as they were for the short nozzle. Note that the transition locations also first occur 
near the corner which is where the above locations are reported. 

Although transition is likely to occur both near the throat and the exit, it is more 
likely to occur near the exit plane where the crossflow Reynolds numbers are very 
large. Therefore, assuming that transition will not occur in the throat, a quiet-test 
Reynolds number can be calculated for the long nozzle. Again, sketching straight lines 
that approximate the characteristics that bound the test rhombus and come from the wall 
approximately 1.0 meter downstream of the throat, a quiet test length of approximately 
0.17m is calculated. This corresponds to a quiet-test Reynolds number (Rax) of 1.2 x 

10^. Although lengthening the nozzle reduces the amount of crossflow apparent in the 
boundary-layer, the relative locations of transition appears to be the same regardless of 
the length. However, these computations of quiet-flow length are very preliminary. 

Although it has been assumed that transition may not occur in the nozzle throat, 
recent computational results [64] suggest that this is not the case. A linear stability 
analysis was performed on the long nozzle geometry using a unit Reynolds number of 
1.0 x 107 ft"l . The results of this analysis state that at 75% of the distance from the 
centerplane to the comer, the N-factor was about 4 or 5 in the throat, 8 at about 0.3 ft, 
and about 9 at 1.0 feet downstream of the throat. In general, the N-factors rise rapidly 
thorough the throat region and then level-off. The results also state that at 89% of the 
way to the corner, the N-factor was 10 at 0.5 ft downstream of the throat. Since N- 
factors greater than 9 usually result in transition, it is very likely that transition will 
occur slightly downstream of the throat in the long nozzle. Although the crossflow 
Reynolds numbers were not accurate in estimating the exact location of transition, the 
data still suggested that transition would occur in the throat region. Therefore, the 
crossflow Reynolds number techniques were helpful in determining the applicability of 
these nozzles for quiet-test nozzles. 


3 .4 Reduction of Crossflow Pressure Gradient in Nozzle Throat 
The two Mach 2.4 nozzle designs provided quiet-test Reynolds numbers that are 
comparable to those in quiet-flow nozzles used today. However, to obtain these results, 
the possibility of transition in the nozzle throat had to be neglected. The linear stability 
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analysis by Malik showed that transition will most likely occur in the throat of the long 
nozzle. Since the crossflow results in the short nozzle are similiar to those in the long 
nozzle, it is reasonable to assume that transition will occur in this nozzle as well. 
Therefore, to design a square nozzle which will not experience transition near the 
nozzle throat, a method must be developed to minimize the crossflow in this region. 
This is the purpose of this section. 

Before new designs can be developed with smaller crossflow in the throat 
region, the reasons for the crossflow must be known. As discussed in the introduction, 
crossflow is created by a pressure gradient that is perpendicular to the freestream 
direction. Therefore, to reduce the crossflow in the nozzle, the pressure gradient 
between the nozzle centerplane and comer must be reduced. In the inviscid design, this 
can be done through the relation 


_P_ 

Pn 


1 + 


I^m 2 1 y-1 

2 J 


where 


y Ratio of Specific Heats 

P Pressure 

P 0 Total Pressure 

M Mach Number 


(37) 


Therefore, the crossflow pressure gradient can be determined as long as the Mach 
numbers are known at the centerplane and the comer. 

Figure 62 shows the Mach number contours at the throat for the Navier-Stokes 
solution in the short nozzle design. As can be seen, the Mach number increases radially 
towards the comer, with the Mach number at the wall in the centerplane being smaller 
than in the comer. Also note that the Mach number at the centerplane and the comer is 
slightly greater than 1.0. This data can then be taken into account by setting the Mach 
number at the comer to be 1.0+e, where e is small. This can then be substituted into the 
above equation, and an approximate solution is found assuming that all terms and 
smaller are neglected. The solution is then found to be 
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Figure 62 -Mach Contours in the Short Nozzle Throat 
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This technique can also be applied in the centerplane of the nozzle by letting M = 1.0 + 
(3 where (3 is small and e > |3. The solution in the centerplane is then found to be 




(39) 


Subtracting these two equations gives the approximate pressure gradient between the 
comer and center of the nozzle 
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where 

6 Mcomer " find is « 1.0 

P Mcenter - 10 and is « 1.0 

e> P 


Now, to make this pressure gradient as small as possible, e - P should go to zero. This 
can be accomplished by making the Mach contours as close to the one-dimensional 
approximation as possible. This means that in the nozzle throat, the sonic line must be 
as vertical as possible. 

This can be analyzed using the relations of Hopkins and Hill [50] that are used 
for the transonic portion of the nozzle design. The equations used in this portion of the 
analysis can be found in detail in reference 52. Referring to figure 63 (figure 1 in 
Hopkins and Hill) the sonic line begins on the axis where x = 0.0. It is also stated that 
on the reference line x = ^, and T) = 0. Therefore, the location of the sonic point is 
coincident with the location where x = ^ =0. Substituting this into H r and its 
derivatives and M r and its derivatives, shown in reference 52, give the following 
relations on the reference line 


H r =1 


H' r =0 
h? = — 

M r = 1 

[ 

Rs 

m H' 

X 

II 

o 

V 

Hf S- y— 

M? = - 


R s C l 


y + 1 

2R7 

(y- 3 Xv-l) 1 
2(y + l) R s 


(41) 


where 


R s = (0.9322 + 0.0565y) 


(0.48Cr° 189 ) 

0.6173e' 1 f + 


R a (l.l234 + 0.0077 IQ - 


(42) 


0.000163CJ) - rJ| (o.0182 + 0.0011 IQ -0. 000020 icj)j 
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Figure 63-Hopltins and Hill [50] Geometry 
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(43) 


Ra Radius of Curvature 

0a Approach Angle 

y Ratio of Specific Heats 


These relations can then be substituted into equation (17d) in Hopkins and Hill to get 
the ratio of dimensionless velocity inside the domain to the dimensionless velocity on 
the reference line, 







(44) 
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This equation can then be evaluated at two locations, q = 1 .0 corresponding to the 
square corner, and q = 1/V2 corresponding to the centerplane. Note that on the 
reference line, M r = Mr = 1.0, simplifying the above equation. The final equations at 
each location are shown below 


M 


= 1 + 

comer 2R 


1 1 ^ 
+ — *-0, 


2Ri 


(45) 


and 


where 


M 


'center 


i 1 1 * 

= 1 + — + — *-<!>, 


4Re 


8R_ 


(46) 


0 = 1 


1- 


2(y + l) 


(47) 


These two equations can be related to the Mach number using the equation 


M 2 = 


2M 


y + 1 - (y - l)M 


„2V 


(48) 


Substitution of equations 45 and 46 into equation 48 gives the Mach number at the 
comer and the wall in the centerplane in terms of the design parameters Ra and 0g. 
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These two equations can then be solved for e and P and substituted back into equation 
37 to give the crossflow pressure gradient in terms of the design variables R a and 0 a . 
This final equation is shown below 


'comer 


-f .r^rAK^w-oi 

0 center L J \Y V 


where 



For (P/Po)comer - (P/Po) center in equation 51 to go to zero, the variable Rs must go to 
infinity. Referring to equation 42, R$ can only go to infinity if R a goes to infinity or if 
Cl goes to zero or infinity. However, equation 43 shows that Cl is bounded to be a 
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finite value by the inverse tangent function. Therefore, the pressure gradient only goes 
to zero if the radius of curvature goes to infinity or the variable C i goes to zero. 

By assuming that y = 1.4, the above equation can be solved for many different 
non-dimensional radii of curvature and approach angles. The results are shown in 
figure 64. For a radius of curvature ratio of less than 2, small changes in the radius of 
curvature have a large effect on the crossflow pressure gradient. After this point, 
however, the change is much more gradual. It should also be noted that the maximum 
radius of curvature used in this analysis is 11.0. After this point, the equations for R s 
and C i break down since they are curve fits initially developed using radii of curvature 
less than 11. 

Figure 65 shows the radii of curvature from 5 to 10. From this it can be seen 
that doubling the radius of curvature gives a savings of less than half in the pressure 
gradient. To verify equation 51, the pressure gradient is estimated from figure 65 at the 
short and long nozzle radii of curvature. For a radius of curvature of 5, like that in the 
short nozzle, the pressure gradient is estimated to be approximately -0.036. This is the 
same result at the throat in figure 25 given by the SNODEC code. For a radius of 
curvature of 6.6 as in the long nozzle, this approximation gives a pressure differential of 
-0.029, the same results as seen in figure 52 given by the SNODEC code. From these 
two results, it seems reasonable to assume that the analysis presented is correct. 

The results from the Navier-Stokes solution, however, are not the same as the 
inviscid design. This could be due to the comer flow or other viscous effects not taken 
into consideration in the inviscid design. Figure 25 shows that the viscous results have 
a smaller pressure differential in the throat region. In fact, the pressure differential is 
approximately 11% smaller than the inviscid design predicted. Figure 52 shows that the 
viscous results also have a smaller differential in the long nozzle. However, the 
difference is slightly less than 9%. Since the difference in pressure gradient declines as 
the radius of curvature increases, it would seem reasonable to assume that as the radius 
of curvature becomes larger, the inviscid approximation becomes closer to the viscous 
results. However, two points of data are not enough to verify this assumption. Further 
analysis is necessary to determine the relation between the crossflow pressure gradient 
in the inviscid design and the viscous results. 

Finally, analysis of how the crossflow pressure gradient affects the crossflow 
itself would be useful for the nozzle design. This cannot be done explicitly, however, 
since no correlations exist relating crossflow to the crossflow pressure gradient. The 
viscous boundary-layer growth is much harder to estimate. The results that have been 
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obtained in this investigation, however, can be used to try to get an idea of how these 
two properties interact. Figures 30 and 32-34 show the WMAXAJe contours and the 
crossflow Reynolds numbers for the short nozzle. Figure 30 shows that the maximum 
value of WMAXAJe the throat is almost 5 corresponding to R c f of a little more than 
500, RcfR&H of approximately 500, and R 52 of about 300. These values occur with a 
crossflow pressure gradient of about -0.0325 as a fraction of the total pressure. In the 
long nozzle, the crossflow pressure gradient is about -0.0265 as a fraction of the total 
pressure, or about 18.5% smaller than in the short nozzle. Figures 53 and 57-59 show 
the WMAXAJe contours and the crossflow Reynolds numbers for the long nozzle. It 
can be seen in figure 53 that the maximum value of WMAXAJe is about 4. This is 20% 

smaller than in the short nozzle. However, in figures 57-59, the maximum crossflow 
Reynolds numbers are R c f of approximately 500, R c fR&H of about 450, and R 52 close 
to 300. This corresponds to only a slight decrease in R c f, a 10% reduction in RcfR&H, 
and no real difference in the value of R 52 • 

This suggests that the decrease in the crossflow pressure gradient has a linear 
effect on the decrease in WmAX- However, with this decrease comes an increase in the 
boundary-layer thickness evident in figure 40. This is most likely the reason for the 
standard crossflow Reynolds number being approximately constant. Therefore, if a new 
nozzle was designed with a radius of curvature of 10 , the crossflow pressure gradient 
would be in the range of -0.021 or about 40% smaller than in the short nozzle design. 
The corresponding maximum value of WMAXAJe in the throat would then be around 
2 %. Since this analysis is very preliminary, it is unclear whether this reduction would 
be enough to prevent transition from occurring in the throat. Further analysis must be 
completed to verify these observations. The effect of crossflow pressure gradient on the 
resulting crossflow boundary-layer profiles is as yet poorly understood. 
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CHAPTER IV CONCLUSIONS 


CFD simulations of the flow through two Mach 2.4 nozzles with square test 
sections have been completed. The initial design of a Mach 8 nozzle with a square test 
section has also been completed. Grid refinement studies have been completed to verify 
the accuracy of the solutions. The results showed that the solutions are accurate enough 
for preliminary design work, but more accurate results must be obtained for reliable 
stability-based transition estimation techniques. The boundary-layer edge defmition 
was also shown to be an important parameter for the calculation of the preliminary 
methods used to estimate transition. 

Three preliminary methods of estimating crossflow-induced transition have been 
applied to the simulation results for both Mach 2.4 nozzles. Estimates have shown that 
transition is likely in the throat, and even more likely near the exit in both nozzles. 
Assuming transition does not occur in the throat of either nozzle, estimates of the quiet 
test length were completed for both nozzles. Results suggest that these nozzles might 
be suitable for quiet-flow facilities. However, recent investigations by other workers 
using e^ techniques suggest that transition occurs slightly downstream of the nozzle 
throat. 

A design procedure has been developed to study the crossflow pressure gradient 
in the nozzle throat. Initial estimates, using the solutions for the Mach 2.4 nozzles, 
suggest that the ratio of radius of curvature to the throat radius of 10 or greater may be 
necessary to eliminate transition in the nozzle throat. These results are preliminary, 
however, and more cases need to be tested to verify this procedure. 
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APPENDIX 


A ppendix A Grid Generation Programs and Data 
The following program is used to convert the *.csd files from the SNODEC 
code to be used in the grid generation program. 

* This program was written to convert data from the *.csd file 

* created by a program called SNODEC for the development of 

* square nozzles. The program reads the data and writes a file 

* for use in a grid generation code. 

* 

* Written by: Tim Alcenius 

* Purdue University 

* June 1993. 

* 

* Identify variabels and open files to be read. 

* 

real x,y,z 
integer i 

open(l, file = 'm8.csd') 
open(2, file = ’pararri) 
xmax = 0. 
xmin = 0. 

* 

* Read the number of data points in any plane. If it is 1, then go to 

* next section of the program. 

* 

20 read(l,*) imax 
if(imax .EQ. 1) then 
goto 30 
endif 

* 

* Read in variables. Set minimum and maximum values of x, and 

* then scale data to metric units. 

* 

do 10 i = 1, imax 
read(l *) y,z,x 
if(x .LT. xmin) then 
xmin = x 
endif 

if(x .GT. xmax) then 
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xmax = x 
endif 

if(z .GT. 0 .AND. y .GT. 0 .AND. x .EQ. xmax) then 
yscale = .3048 
endif 

10 continue 
goto 20 

* 

* Rewind the file for which data was read and write the mimimum and 

* maximum values of x to the output file. 

* 

30 rewind(l) 

write(2,*) xmin*yscale, yscale*xmax 
print *, yscale 

* 

* Read the number of points in the current data set. If the value is 

* 1, then quit. 

* 

40 read(l,*) imax 
if(imax .EQ. 1) then 
goto 50 
endif 

* 

* Write data to the output file in the form of x, ymin, ymax, zmin, zmax. 

* 

do 60 i = 1, imax 
read(l,*) y,z,x 

if(y .EQ! z .AND. y .GT. 0 .AND. z .GT. 0 .AND. i .LT. imax) then 
write(2,*) x*yscale,a,y*yscale,b,z*yscale 
endif 

60 continue 
goto 40 
50 end 


The following file is the *.csd file for the short Mach 2.4 nozzle created by the 
SNODEC code. This is given as an example of this type of file. 


0.253884 

0.253884 

-0.253884 

-0.253884 

0.253884 

5 

0.251907 

0.251907 

-0.251907 


0.253884 

-0.253884 

-0.253884 

0.253884 

0.253884 

0.251907 

-0.251907 

-0.251907 


-0.437861 

-0.437861 

-0.437861 

-0.437861 

-0.437861 

-0.426026 

-0.426026 

-0.426026 
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0.251907 -0.426026 
0.251907 -0.426026 


-0.251907 

0.251907 

5 

0.249963 

0.249963 

-0.249963 

-0.249963 

0.249963 


0.334358 

0.334358 

-0.334358 

-0.334358 

0.334358 

5 

0.334391 

0.334391 

-0.334391 

-0.334391 

0.334391 

5 

0.334421 

0.334421 

-0.334421 

-0.334421 

0.334421 

1 


0.249963 

-0.249963 

-0.249963 

0.249963 

0.249963 


0.334358 

-0.334358 

-0.334358 

0.334358 

0.334358 

0.334391 

-0.334391 

-0.334391 

0.334391 

0.334391 

0.334421 

-0.334421 

-0.334421 

0.334421 

0.334421 


-0.414193 

-0.414193 

-0.414193 

-0.414193 

-0.414193 


2.35498 

2.35498 

2.35498 

2.35498 

2.35498 

2.36681 

2.36681 

2.36681 

2.36681 

2.36681 

2.37865 

2.37865 

2.37865 

2.37865 

2.37865 


The following two programs are used in the initial generation of the nozzles. 
The first is called sivinput and is used to generate the data used in the Sivells code. The 
second is entitled makeblin and is used to make the input file for the Harris boundary - 
layer code. 

* From moen Wed Nov 25 12:00:09 1992 

c this program writes an input file in the correct format for 
c sivells code; sps 6-27-90 

c modified 9-30i-92 for the square nozzle work with streamlines sps 
c modified to produce starting streamlines on nozzle wall. 

* mod sps 1 1-30-92 to make mp=5 

* mod 1 1-30 sps to allow choices of IN, XC, IX 
c 

character* 10 title 
character*20 sivfile 
c 
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write (*,*) 'enter a rootname to write sivells input file:' 
read (*,10) sivfile 
title = sivfile 

ileng = index(sivfile,' ')-l 
sivfile(ileng+l:ileng+4) = '.inp' 
write (*,*) 'opening file-',sivfHe,'-for output' 
open(unit=2,file=sivfile, status- new') 
c 

* write (*,*) 'enter title of run (10 characters): ' 

* read (*,20) title 

write (*,*) 'enter jd (-1=2D, 0=axisym): ' 
read (*,*) jd 
write (2,30) tide,jd 

* write (*,*) 'enter sfoa:' 

* read (*,*) sfoa 

sfoa=0. !use 3rd or 4th degree distribution 
gam = 1 .40 
ar= 1716.563 
zo=l 

c following three used in bl computations, not used here 
ro=l 
visc=l 
vism=l 
c 

xbl=1000. ! gives values at evenly spaced intervals 
write (2,40) gam,ar,zo/o,visc,vism,sfoa,xbl 
c 

write (*,*) 'enter etad,rc,bmach,cmc: ' 

read (*,*) etad,rc,bmach,cmc 

xc=0. Iso 4th degree distribution, change? 

* write (*,*) 'enter xc, in: ' 

* read (*,*) xc,in 

* write (*,*) ’xc,in= ',xc,in 

fmach=0. I this sets distribution, change? 
sf = 0. Inozzle throat radius = 1.0 
pp = 0. [coordinates given relative to throat 
write (2,50) etad/c,fmach,bmach,cmc,sf,pp,xc 
c write (*,*) 'enter mt,nt,ix,in,md,nd,nf^np,jc,lr,nx:' 
c read (*,*) mt,nt,ix,in,md,nd,nf,mp,jc,lr,nx 
mt=61 Ipts on char. EG, max 125 
nt=31 Ipts on axis IE, max 149-LR 

* write (*,*) 'enter ix: ' 

* read (*,*) ix 

ix=0 I is 3rd deriv matched? change? 

in=10 fuse Mach no. distrib on BC, makes 2nd deriv match rad. flow 
c change? 

iq=0 I calls for complete contour 
md=61 Ipts on char. AB, max about 125, odd 
nd=15 Ipts on axis BC, max about 150, 
c changed from 49 to 15 sps 7-2-90 

write (*,*) 'enter -1 for smoothing, 1 for no smoothing: ' 
read (*,*) ismooth 
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nf=ismooth*8 1 !pts on characteristic CD. Neg calls for smoothing 

mp=5 (pts on GA, conical section, if Fmach ne Bmach 

jc=0 (if not 0, used to print intermediate characteristics 

lr=31 !pts on throat char., - prints out transonic soln 

nx=13 (spacing of pts on axis upstream, this no. recc. 

mq=0 (pts downstream of D 

jb=- 1 ! neg for no BL computation 

jx=l (pos calls for streamlines 

it=0 (jack points, not used 

write (2,60) mt,nt,ix,in,iq,md,nd,nf^np/nqjb,jx,jc, 

> it,lr,nx 

if (ismooth .eq. -1) then 

noup= 1 0 ! smoothing parameters, arbitrary 

nodo=10 

npct=90 

write (2,70) noup,npct,nodo 
end if 


c 

c gives streamline distribution that corresponds 
c to the half wall for conversion to a square nozzle, 
c note that the number of streamlines requested will 
c be reduced by one because Sivells automatically 
c calculates the wall streamline. Sivells output 
c will have the actual number of streamlines requested, 
c (moen 10-92) 
c 

write (*,*) 'How many streamlines along halfwall?' 

read (*,*) nstream 

nstream=nstream- 1 

dx= 1 .0/(float(nstream)*sqrt(2.0)) 

ycnt=1.0/sqrt(2.0) 

do 100 istream = Onstream- 1 

etadstr = etad*sqrt((istream*dx)**2+ycnt**2) (see btm page 59 
qm = sqrt((istream*dx)**2+ycnt**2) 
xj = 1 (look for more streamlines 
write (2,90) etadstr ,qm,xj 
if (ismooth .eq. -1) then 

write (2,70) noup,npct,nodo (must have for each! ! 
end if 

100 continue 


close(2) 

stop 

10 format(a20) 

20 format(alO) 

30 format(lx,al0,2x,i2) 
40 format(8(lx,f9.3)) 

50 format(8(lx,f9.3)) 

60 format( 1 x,i4, 1 5(i5)) 
70 format(lx44,2(i5)) 

90 format(3(lx,f9.4),lx) 
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end 

* SUN VERSION WITH SPACES BEFORE ITEMS IN NAMELIST ! ! 9-92 

* Steven P. Schneider Purdue University 3 17-494-3343 

* this is a program to read in output from the sivells code 

* add specifics for Re, and write in a form readable 

* by the Harris code for bl. 

* specific for the nozzle block problem sps 6-90 

* add some code for the contraction computation 12-5-90 sps 

* add code for output of file for arbitrary shape using modified Newtionian thy 

* sps 3-6-91 

* allow for arbitrary power-law distribution of points 3-6-91 sps, 

* and for easy change of jsolve 

* add code for reading in pressure data from euler code, and shock location 

* data, and writing an output file for the harris code 7-13-91 sps 

* modified to fix bug with error code on different case 7-28-91 sps 

* 

parameter(maxpts= 1000,jsolve= 1 ) 

* jsolve is the number of solution stations per station of input data 

* can increase to make finer resolution solution possible 
character*80 text 

character* 20 root,infile,outfile,parmfl 

* most of these are the number of stations along the moc soln 

* ss is the number of stations in bl soln to get 

dimension x(maxpts),y(maxpts),pepo(maxpts),s(maxpts),theta(maxpts) 
dimension proval(maxpts),pmtval(maxpts),ss(jsolve*maxpts) 
dimension tw(maxpts) 

dimension xsh(maxpts),ysh(maxpts) Ishock location 

common /param/pstar !for passing to subroutines for arbitrary shapes 

data pi/3.1415926535/,kspmt/l/,ksprof/5/ 

* * (print info every kspmt’th soln station; print profile info 

* * (every ksprof th soln station) 

* need dense printing of soln for gortler work to get good values 

* for derivatives of wall height to get streamwise curvature! 

write (*,*) 'enter root filename for transfer: ’ 
read (*,5) root 
5 format(a20) 
ileng = index(root,' ') - 1 
write (*,*) 'enter 0 if this is a Sivells test, ' 
write (*,*) 'enter 1 if this is a flat plate test, ' 
write (*,*) 'enter 2 if this is a Lees modified newtonian test: ' 
write (*,*) 'enter 3 if Euler output for body is to be read: ' 
read (*,*) imodel 
if (imodel .eq. 0) then 
write (*,*) 'Sivells test' 
infile(l:ileng) = root(l: ileng) 
infile(ileng+l:ileng+3) = '.bl' 
write (*,*) 'reading input data from file '.infile 
open(unit= 1 ,file=infile, status- old') 
read (1,10) text 
10 format(a80) 



read (1,10) text 
read (1,*) numpts 

if (numpts .gt. maxpts) stop 'too many points' 
do 50 i=l, numpts 
read (1,*) x(i),y(i),pepo(i) 

50 continue 
close (uni t=l) 

(done reading in info from sivells output file) 
else if (imodel .eq. 1) then 
write (*,*) 'enter numpts, mach, gam, plend: ' 
read (*,*) numpts, amach,gam,plend 
gmexp = gam/(gam- 1) 
gmfact = (gam-l)/2.0 

denom = (1.0 + gmfact*amach**2)**gmexp 
pepO = 1.0/denom 
write (*,*) ' gives pep0= '.pepO 
do 60 i = 1, numpts 

x(i) = plend*float(i-l)/float(numpts-l) 
y(i) =1.0 !not 0.0, messes up computations 
pepo(i) = pepO 
60 continue 

else if (imodel .eq. 2) then 
write (*,*) 'this is a modified newtonian test;' 
write (*,*) 'you must enter the shape in source code' 
write (*,*) 'enter numpts, mach, pstar, plend: ' 
read (*,*) numpts, amach,pstar,plend 
gam =1.4 lair 

if (numpts .gt. maxpts) stop 'too many points' 
now compute pressure ahead of shock, ratio to total pressure 
in stilling chamber 
gmexp = gam/(gam-l) 
gmfact = (gam-l)/2.0 

denom = (1.0 + gmfact*amach**2)**gmexp 
pinfpO = 1.0/denom 
write (*,*) ' gives pinfinity/p0= ', pinfpO 
now compute stagnation or total pressure behind normal shock, 
ratio to p_infty ahead of shock(see Anderson p. 54, 3.17) 
denoml = 4*gam*amach**2 - 2*(gam-l) 
piecel = ( (gam+l)**2 * amach**2/denoml )**(gmexp) 
piece2 = ( 1 -gam+2 *gam*amach* *2)/(gam+ 1 ) 
ptpinf = piece l*piece2 
ptpO = ptpinf*pinfpO 

write (*,*) 'pstag on body/p infty is '.ptpinf 
write (*,*) 'gives pstag on body/pstag in stilling=',ptp0 
cpmax = (2.0/(gam*amach**2)) * (ptpinf - 1) ! Anderson 3.19 
write (*,*) 'cpmax computed as '.cpmax 
now compute body shape 
for non-sphere shape, do this in subroutine 
write (*,*) 'passing to blunts; pstar ,plend=',pstar,plend 
pause 'need to incorporate subroutine with body shape! !' 
call blunts(plend, numpts, x,y, theta) 
do 80 i = 1, numpts 


* now compute ratio of pe to ptotal ahead of shock, using Lees 

* modified newtonian thy - formula derived using p. 54 

pepo(i) = (ptpO - pin^>0)*(sin(theta(i)))**2 + pinfpO 
80 continue 
else if (imodel .eq. 3) then 
write (*,*) 'working for euler data for blunt body' 
write (*,*) 'enter mach, gamma: ' 
read (*,*) amach,gam 

* now compute pressure ahead of shock, ratio to total pressure 

* in stilling chamber 

gmexp = gam/(gam-l) 
gmfact = (gam-l)/2.0 

denom = (1.0 + gmfact*amach**2)**gmexp 

pinfpO = 1 .0/denom 

write (*,*) ' gives pinfinity/pO= '.pinfpO 

* now compute stagnation or total pressure behind normal shock, 

* ratio to p_infty ahead of shock(see Anderson p. 54, 3.17) 

denoml = 4*gam*amach**2 - 2*(gam-l) 

piecel = ( (gam+l)**2 * amach**2/denoml )**(gmexp) 

piece2 = ( 1 -gam+2 *gam* amach* * 2)/(gam+ 1 ) 

ptpinf = piece l*piece2 

ptpO = ptpinf*pinfpO 

write (*,*) 'pstag on body/p infty is '.ptpinf 

write (*,*) 'gives pstag on body/pstag in stilling=',ptpO 

write (*,*) 'reading euler data for blunt body: ' 

infile( 1 :ileng) = root( 1 :ileng) 

infile(ileng+l:ileng44) = ’.eul' 

write (*,*) 'reading input data from file '.infile 

open(unit= 1 ,file=mfile,status- old') 

read (1,90) lineskip 

90 format(i6) [number of lines to skip 

do 92 i = 1, lineskip [skip 'lineskip' lines of text 
read (1,91) text 

91 format(a80) 

92 continue 

read (1,*) numpts 

if (numpts .gt. maxpts) stop 'too many points' 
do 95 i=l, numpts 

read (1,*) x(i),y(i),pepinf [read euler data for pressure 
pepo(i) = pepinf*pinfpO 

95 continue 

* Now read shock location data: 

read (1,90) lineskip 

do 96 i — 1, lineskip [skip 'lineskip' lines of text 
read (1,91) text 

96 continue 

read (1,*) numshpts 

if (numshpts .gt. maxpts) stop 'too many points' 
do 97 i=l, numshpts 

read (1,*) xsh(i),ysh(i) [read euler data for shock location 

97 continue 

write (*,*) 'done reading from file’ 



else 

stop 'imodel must be 0,1,2, or 3' 
end if 

* now check computed values 

do 150 i = l.numpts 
if (imodel .gt. 1) then 

if (pepo(i)-ptpO .gt. -0.01*ptp0) then ! arbitrary nearness 
write (*,*) 'pepo(',i,')= ’,pepo(i) 
write (*,*) 'ptpO= ',ptpO 

write (*,*) 'if first gt second when computed by VGBLP,' 
write (*,*) 'this will give surface static pressure larger' 
write (*,*) ’ than the total pressure on the surface' 
write (*,*) 'this would be a fatal error in VGBLF 
if (pepo(i) .gt. ptpO) then 
write (*,*) 'reducing pepo(’ 4 ,') to 0.99999*ptp0' 
write (*,*) ' to forstall error in VGBLP' 
pepo(i) = 0.99999*ptp0 
end if 
end if 
end if 

150 continue 

* 

* now read in parametric info from parameter file 

* 

parmfl(l:ileng) = root(l:ileng) 
parmfl(ileng+l:ileng+3) = '.re' 

write (*,*) 'reading reynolds number scaling info from '.paimfl 

open(unit= 1 ,file=parmfl,status- old') 

read (1,5) outfile !read filename to write to 

read ( 1 ,*) throat ! throat radius in feet 

* (assumes input scaled so throat radius is one unit) * 

read (1,*) prandtl 

* ptotal,ttotal,xmach are conditions at infinity - see p. 44 

* for nozzle are stagnation chamber conditions 

* for non-nozzle, xmach seems to affect mostly the computations 

* involving the flow behind the shock. Should be freestream values! 
read (1,*) ptotal 

read (1,*) rgas 

read (1 ,*) ttotal ! ahead of le shock 
read (1,*) xmachi 
if (imodel .ne. 0) then 
if (xmachi .It. 1.0) then 
write (*,*) 'xmach given as '.xmachi 
write (*,*) 'should be freestream value ahead of shock,' 
write (*,*) ' not the value at stagnation!!' 
end if 
end if 

close (unit=l) 

* 

* now set defaults for input to harris program (besides harris's) 

* 

if (imodel .eq. 2 .or. imodel .eq. 3) then 
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ibody = 1 !stag pt at nose 
j = 1 laxisymmetric 

else 

ibody = 2 !no stagnation point at nose 
j = 0 !2D 

end if 

ie = 5 1 Iffom test case number 4 
if (imodel .eq. 3) then 
ientro = 2 ! variable entropy calculation 
else 

ientro = 1 
end if 

igeom = 1 ! create coords using geometric series; need xk,ie,xend! ! 

kodunit=0 !US units 
if (imodel .eq. 0) then 

kodwal=l ! specify wall temperature distribution (no time to heat) 
else 

kodwal=2 ! specify adiabatic wall 
end if 

iendl = jsolve*numpts Inumber of soln stations 
proinc = 10.0 ! hopefully, none 
pmtinc = 10.0 

sst = le+20 !no transition on body until then (laminar flow) 
if (imodel .gt. 1) then 

wave = 90.0 ! shock wave angle at s=0, needed for test case 4 type flows 
else 

wave = 0.0 (needed for shockless type flows 
end if 

c xend = 120 Iffom blasius test case 
xend =10 las in test case 4 

c xk = 1.275 lvalue used in test cases in book, sets grid 
c xk = 1 . 1 (because value used in text gives hyper-dense 

* grid near wall, which makes for difficulties, 
c xk = 1.0 Hike test case number 4 

c xk = 1.05 I because 1.0 gives little near wall for sphere 

xk = 1 . 1 (because 1 .05 gives not great resolution for stability 

* 

* compute the arc length along the wall (see (66) of harris paper) 

* (approximate with straight line segments between stations) 
s(l) = 0.0 

write (*,*) 'setting wall temp equal to total temp, beware!' 

* model temp same as total temp, ambient 

* note that total temperature behind shock is the same as ahead of shock 

* *(true at stagnation point, anyway)* 
tw(l) = ttotal 

do 100 i = 2,numpts 

s(i) = s(i-l) + sqrt( (x(i)-x(i-l))**2 + (y(i)-y(i-l))**2 ) 
tw(i) = ttotal I set above, assumes same along wall 
100 continue 

* compute the s stations to get soln at, and to write at, even spacing 

* in sqrt(s) normally, other power for other cases 
if (imodel .It. 2) then 


power = 2.0 

wnte (*,*) 'using square root distribution of pts' 
else 

* with sphere, have problems with stepsize increasing too rapidly near le 

power = 1.0 

write (*,*) 'using linear distribution of pts, good for blunt' 

* *(try this, since problems with T<0 at le)*like test case 4 
end if 

rootl = 1.0/power 

rsinc = (s(numpts)-s( 1 ))* *root l/float(iend 1 ) (try- 

ss(l) = rsinc* *power 

ss(2) = ss(l) (required 

ss(3) = ss(l) (also required, actually 

srun = 3*ss(l) (running value of s 

ipmt = 0 

ipro = 0 

do 200 i = 4,iendl 
srunold = srun 

srun = ( (srunold) **rootl + rsinc )**power 
ss(i) = srun-srunold 

if (mod(i,kspmt) .eq. 0 .or. i .eq. 1) then 
ipmt = ipmt + 1 
pmtval(ipmt) = srun 
end if 

if (mod(i,ksprof) .eq. 0) then 
ipro = ipro + 1 
proval(ipro) = srun 
end if 

200 continue 

* 

* Now write out NAM1 namelist into file 

* 

write (*,*) ' writing bl input to file '.outfile 

open(unit=2,file=outfile) 

write (2,*) ' &NAM1' 

write (2,*) ' IBODY-.ibody 

write (2,*) ' IE=',ie 

write (2,*) ’ ENDl=',iendl 

wnte (*,*) 'may need to reset param.inc, note that ' 

write (*,*) 'ie is now ',ie,' and iendl= ',iendl 

write (2,*) ’ IENTRO=',ientro 

write (2,*) ' IGEOM=',igeom 

write (2,*) ’ EPRO-.ipro 

write (2,*) ' IPRNT-.ipmt 

write (*,*) 'ipro and ipmt are ',ipro,ipmt 

write (2,*) ’ J=’,j 

write (2,*) ' KODUNIT-.kodunit 
write (2,*) ' KODWAL=',kodwal 
if (imodel .ne. 2) then 

phii = (180.0/pi)*atan( (y(2)-y(l))/(x(2)-x(l» ) 
else 

write (*,*) 'setting leading edge angle = 90 degrees’ 
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phii = 90.0 
end if 

write (2,*) ’ PHII=',phii 

write (2,*) ' PR=',prandtl 

write (2,*) ' PRNTINC-.pmtinc 

write (2,*) ' PRNTVAL=',(prntval(i)*throat,i=l,iprnt) 

write (2,*) ' PROINC- .proinc 

write (2,*) ' PRO V AL=' ,(proval(i) * throat, i= 1 ,ipro) 

write (2,*) ' PTl=',ptotal 

write (2,*) ' R=',rgas 

write (2,*) ' SST=',sst 

write (2,*) ' TTl=',ttotal 

write (2,*) ' WAVE=',wave 

write (2,*) ' XEND=’,xend 

write (2,*) ' XK-,xk 

write (2,*) ' XMA=',xmachi 

write (2,*) ' &END' lend of namelist input 

* end of writing naml. Now compute and write nam2: 

write (2,*) ' &NAM2' 
write (2,*) ’ NUMBER- .numpts 
write (*,*) 'and NUMBER is '.numpts 
write (2,*) ' PE=',(pepo(i)*ptotal,i=l,numpts) 
write (2,*) ' RMI=' ,(y (i)* throat,i= 1 .numpts) 
write (2,*) ' S- ,(s(i)*throat,i=l, numpts) 
write (2,*) ' SS=',(ss(i)*throat,i=l,iendl) 
if (imodel .eq. 0) then 
write (2,*) ' TW- ,(tw(i),i=l,numpts) 
else 

write (2,500) numpts- 1 'make adiabatic wall 

500 format(' QW- ,i4,'*0.0,0.0') 1,0.0 kluge few compiler bug 
end if 

write (2,*) ' Z- ,(x(i)*throat,i=l,numpts) 
write (2,*) ' &END' 

* end of writing nam2. Now write nam3, if required 

if (imodel .eq. 3) then ! doing entropy computations through shock 
write (2,*) '&NAM3' 
write (2,*) 'NUMBER-, numshpts 
write (*,*) 'now numshpts (JL) is '.numshpts 
write (2,*) 'RRS= ',(xsh(i)*throat,i=l, numshpts) 
write (2,*) 'ZZS= ’,(ysh(i)*throatJi=l,numshpts) 
end if 

write (2,*) 7 
close (unit=2) 

* 

stop 'end of makeblin' 
end 


The following file is created by the program described above. The short Mach 
2.4 nozzle data is shown as an example of this type of file. This file is labeled PAR AM 
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and is the exact coordinates for the short Mach 2.4 nozzle. The first line is minimum 
and maximum x values. Then the values are x, minimum y, maximum y, minimum z, 
and maximum z for one quadrant of the nozzle. Therefore, at any x location, the 
contour in the centerplane is found by taking the minimum z value and the maximum y 
value. 


-0.133460 0.725013 

-0.133460 0. 7.73838E-02 0. 7.73838E-02 
-0.129853 0. 7.678 13E-02 0. 7.678 13E-02 
-0.126246 0. 7.61887E-02 0. 7.61887E-02 
-0.122639 0. 7.56108E-02 0. 7.56108E-02 
-0.119032 0. 7.50485E-02 0. 7.50485E-02 
-0.115425 0. 7.45017E-02 0. 7.45017E-02 
-0.111818 0. 7.39707E-02 0. 7.39707E-02 
-0.108211 0. 7.34556E-02 0. 7.34556E-02 
-1.04604E-01 0. 7.29563E-02 0. 7.29563E-02 
-1.00997E-01 0. 7.24735E-02 0. 7.24735E-02 
-9.73897E-02 0. 7.20066E-02 0. 7.20066E-02 
-9.37827E-02 0. 7.15555E-02 0. 7.15555E-02 
-9.01757E-02 0. 7.11208E-02 0. 7.11208E-02 
-8.65687E-02 0. 7.07023E-02 0. 7.07023E-02 
-8.29614E-02 0. 7.03000E-02 0. 7.03000E-02 
-7.93544E-02 0. 6.99138E-02 0. 6.99138E-02 
-7.57474E-02 0. 6.95441E-02 0. 6.95441E-02 
-7.21404E-02 0. 6.91905E-02 0. 6.91905E-02 
-6.85334E-02 0. 6.8853 IE-02 0. 6.8853 IE-02 
-6.49264E-02 0. 6.85321E-02 0. 6.85321E-02 
-6.13194E-02 0. 6.82277E-02 0. 6.82277E-02 
-5.77124E-02 0. 6.79396E-02 0. 6.79396E-02 
-5.41054E-02 0. 6.76677E-02 0. 6.76677E-02 
-5.04983E-02 0. 6.74126E-02 0. 6.74126E-02 

-4.68913E-02 0. 6.71737E-02 0. 6.71737E-02 
-4.32843E-02 0. 6.695 UE-02 0. 6.695 11E-02 
-3.96773E-02 0. 6.67454E-02 0. 6.67454E-02 
-3.60703E-02 0. 6.65561E-02 0. 6.65561E-02 
-3.24630E-02 0. 6.63833E-02 0. 6.63833E-02 
-2.88561E-02 0. 6.62272E-02 0. 6.62272E-02 
-2.52492E-02 0. 6.60880E-02 0. 6.60880E-02 
-2.16421E-02 0. 6.59651E-02 0. 6.59651E-02 
-1.8035 IE-02 0. 6.58594E-02 0. 6.58594E-02 
-1.44281E-02 0. 6.57701E-02 0. 6.57701E-02 
-1.0821 IE-02 0. 6.56975E-02 0. 6.56975E-02 
-7.21404E-03 0. 6.56417E-02 0. 6.56417E-02 
-3.60700E-03 0. 6.56015E-02 0. 6.56015E-02 
4.41978E-08 0. 6.55820E-02 0. 6.55820E-02 
3.60703E-03 0. 6.55835E-02 0. 6.55835E-02 
7.21407E-03 0. 6.55945E-02 0. 6.55945E-02 
1.0821 IE-02 0. 6.56207E-02 0. 6.56207E-02 
1.44281E-02 0. 6.56640E-02 0. 6.56640E-02 
1.80351E-02 0. 6.57234E-02 0. 6.57234E-02 



2.16422E-02 
2.52492E-02 
2.88562E-02 
3.24633E-02 
3.60703E-02 
3.96773E-02 
4.32843E-02 
4.689 13E-02 
5.04983E-02 
5.41054E-02 
5.77124E-02 
6.13194E-02 
6.49264E-02 
6.85334E-02 
7.21407E-02 
7.57477E-02 
7.93547E-02 
8.296 17E-02 
8.65687E-02 
9.01757E-02 
9.37827E-02 
9.73897E-02 
1.00997E-01 
1.04604E-01 
0.108211 0. 
0.111818 0. 
0.115425 0. 
0.119032 0. 
0.122639 0. 
0.126246 0. 
0.129853 0. 
0.133460 0. 
0.137067 0. 
0.140674 0. 
0.144281 0. 
0.147888 0. 
0.151495 0. 
0.155102 0. 
0.158709 0. 
0.162316 0. 
0.165923 0. 
0.169530 0. 
0.173137 0. 
0.176744 0. 
0.180351 0. 
0.183958 0. 
0.187565 0. 
0.191172 0. 
0.194779 0. 
0.198386 0. 
0.201993 0. 
0.205600 0. 


0. 6.57990E-02 

0. 6.58898E-02 
0. 6.59959E-02 
0. 6.61 163E-02 
0. 6.62504E-02 
0. 6.63973E-02 

0. 6.65564E-02 
0. 6.67268E-02 
0. 6.6907 6E-02 
0. 6.7098 IE-02 
0. 6.72974E-02 
0. 6.75047E-02 
0. 6.77195E-02 
0. 6.79414E-02 
0. 6.81697E-02 
0. 6.84038E-02 

0. 6.86434E-02 
0. 6.88878E-02 
0. 6.91369E-02 
0. 6.93902E-02 
0. 6.9647 IE-02 

0. 6.99074E-02 
0. 7.0171 IE-02 
0. 7.04378E-02 
7.07069E-02 0. 
7.09785E-02 0. 
7.12522E-02 0. 
7.15277E-02 0. 
7.18051E-02 0. 
7.20840E-02 0. 
7.23644E-02 0. 
7.26463E-02 0. 
7.29289E-02 0. 
7.32130E-02 0. 
7.34976E-02 0. 
7.37829E-02 0. 
7.40691E-02 0. 
7.43560E-02 0. 
7.4643 IE-02 0. 
7.49308E-02 0. 
7.52188E-02 0. 
7.55072E-02 0. 
7.57958E-02 0. 
7.60848E-02 0. 
7.63737E-02 0. 
7.66627E-02 0. 
7.69522E-02 0. 
7.72415E-02 0. 
7.7531 IE-02 0. 
7.78203E-02 0. 
7.81102E-02 0. 
7.83994E-02 0. 


0. 6.57990E-02 
0. 6.58898E-02 
0. 6.59959E-02 
0. 6.61163E-02 
0. 6.62504E-02 
0. 6.63973E-02 
0. 6.65564E-02 
0. 6.67268E-02 

0. 6.69076E-02 

0. 6.7098 IE-02 
0. 6.72974E-02 
0. 6.75047E-02 
0. 6.77 195E-02 
0. 6.794 14E-02 
0. 6.81697E-02 
0. 6.84038E-02 
0. 6.86434E-02 
0. 6.88878E-02 
0. 6.91369E-02 
0. 6.93902E-02 
0. 6.9647 IE-02 
0. 6.99074E-02 
0. 7.017 11E-02 
0. 7.04378E-02 
7.07069E-02 
7.09785E-02 
7.12522E-02 
7.15277E-02 
7.18051E-02 
7.20840E-02 
7.23644E-02 
7.26463E-02 
7.29289E-02 
7.32130E-02 
7.34976E-02 
7.37829E-02 
7.40691E-02 
7.43560E-02 
7.4643 IE-02 
7.49308E-02 
7.52188E-02 
7.55072E-02 
7.57958E-02 
7.60848E-02 
7.63737E-02 
7.66627E-02 
7.69522E-02 
7.72415E-02 
7.753 11E-02 
7.78203E-02 
7.81102E-02 
7.83994E-02 
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0.209207 0. 7.86887E-02 0. 7.86887E-02 
0.212814 0. 7.89779E-02 0. 7.89779E-02 
0.216421 0. 7.92672E-02 0. 7.92672E-02 
0.220028 0. 7.95565E-02 0. 7.95565E-02 
0.223635 0. 7.98457E-02 0. 7.98457E-02 
0.227243 0. 8.01347E-02 0. 8.01347E-02 
0.230850 0. 8.04239E-02 0. 8.04239E-02 
0.234457 0. 8.07129E-02 0. 8.07129E-02 
0.238064 0. 8.10021E-02 0. 8.10021E-02 
0.241671 0. 8.1291 IE-02 0. 8.12911E-02 
0.245278 0. 8.15803E-02 0. 8.15803E-02 
0.248885 0. 8.18699E-02 0. 8.18699E-02 
0.252492 0. 8.21598E-02 0. 8.21598E-02 
0.256099 0. 8.24496E-02 0. 8.24496E-02 
0.259706 0. 8.27398E-02 0. 8.27398E-02 
0.263313 0. 8.30297E-02 0. 8.30297E-02 
0.266920 0. 8.33195E-02 0. 8.33195E-02 
0.270527 0. 8.36094E-02 0. 8.36094E-02 
0.274134 0. 8.38986E-02 0. 8.38986E-02 
0.277741 0. 8.41876E-02 0. 8.41876E-02 
0.28134 8 0. 8.44762E-02 0. 8.44762E-02 
0.284955 0. 8.47643E-02 0. 8.47643E-02 
0.288562 0. 8.50520E-02 0. 8.50520E-02 
0.292169 0. 8.53388E-02 0. 8.53388E-02 
0.295776 0. 8.56250E-02 0. 8.56250E-02 
0.299383 0. 8.59103E-02 0. 8.59103E-02 
0.302990 0. 8.61944E-02 0. 8.61944E-02 
0.306598 0. 8.64776E-02 0. 8.64776E-02 
0.310204 0. 8.67592E-02 0. 8.67592E-02 
0.313810 0. 8.70393E-02 0. 8.70393E-02 
0.317419 0. 8.73176E-02 0. 8.73176E-02 
0.321025 0. 8.75940E-02 0. 8.75940E-02 
0.324633 0. 8.78684E-02 0. 8.78684E-02 
0.328239 0. 8.81405E-02 0. 8.81405E-02 
0.331845 0. 8.84106E-02 0. 8.84106E-02 
0.335454 0. 8.86779E-02 0. 8.86779E-02 
0.339060 0. 8.8943 IE-02 0. 8.8943 IE-02 
0.342668 0. 8.92058E-02 0. 8.92058E-02 
0.346274 0. 8.94661E-02 0. 8.94661E-02 
0.349880 0. 8.97237E-02 0. 8.97237E-02 
0.353489 0. 8.99785E-02 0. 8.99785E-02 
0.357095 0. 9.02309E-02 0. 9.02309E-02 
0.360703 0. 9.04805E-02 0. 9.04805E-02 
0.364309 0. 9.07274E-02 0. 9.07274E-02 
0.367918 0. 9.097 15E-02 0. 9.09715E-02 
0.371524 0. 9.12129E-02 0. 9.12129E-02 
0.375130 0. 9.14516E-02 0. 9.14516E-02 
0.378738 0. 9.16872E-02 0. 9.16872E-02 
0.382344 0. 9.19201E-02 0. 9.19201E-02 
0.385953 0. 9.21502E-02 0. 9.21502E-02 
0.389559 0. 9.23776E-02 0. 9.23776E-02 
0.393165 0. 9.26019E-02 0. 9.26019E-02 


0.396773 0 
0.400379 0 
0.403988 0 
0.407594 0 
0.411200 0 
0.414808 0 
0.418414 0 
0.422023 0 
0.425629 0 
0.429235 0 
0.432843 0 
0.436449 0 
0.440058 0 
0.443664 0 
0.447270 0 
0.450878 0 
0.454484 0 
0.458093 0 
0.461699 0 
0.465308 0 
0.468913 0 
0.472519 0 
0.476128 0 
0.479734 0 
0.483343 0 
0.486948 0 
0.490554 0 
0.494163 0 
0.497769 0. 
0.501378 0. 
0.504983 0. 
0.508589 0. 
0.512198 0. 
0.515804 0. 
0.519413 0. 
0.523019 0 
0.526624 0. 
0.530233 0. 
0.533839 0. 
0.537448 0. 
0.541054 0. 
0.544659 0. 
0.548268 0. 
0.551874 0. 
0.555483 0. 
0.559089 0. 
0.562697 0. 
0.566303 0. 
0.569909 0. 
0.573518 0. 
0.577124 0. 
0.580732 0. 


9.28235E-02 0 
9.30420E-02 0 
9.32575E-02 0 
9.34703E-02 0 
9.36800E-02 0 
9.38869E-02 0 
9.40908E-02 0 
9.42917E-02 0 
9.44898E-02 0 
9.46852E-02 0 
9.48775E-02 0 
9.50668E-02 0 
9.5253 IE-02 0 
9.54365E-02 0 
9.56173E-02 0 
9.57947E-02 0 
9.59696E-02 0 
9.614 15E-02 0 
9.63104E-02 0 
9.64765E-02 0 
9.66399E-02 0 
9.68002E-02 0 
9.69575E-02 0 
9.71123E-02 0 
9.72641E-02 0 
9.74132E-02 0 
9.75595E-02 0 
9.77027E-02 0 
9.78436E-02 0 
9.798 13E-02 0. 
9.81163E-02 0. 
9.8248 6E-02 0. 
9.83785E-02 0. 
9.85053E-02 0. 
9.86296E-02 0. 
9.87509E-02 0. 
9.88698E-02 0. 
9.89859E-02 0. 
9.90996E-02 0. 
9.92106E-02 0. 
9.93191E-02 0. 
9.94248E-02 0. 
9.95279E-02 0. 
9.96284E-02 0. 
9.97266E-02 0. 
9.98223E-02 0. 
9.99156E-02 0. 
1.00006E-01 0. 
1.00094E-01 0. 
1.00180E-01 0. 
1.00264E-01 0. 
1.00344E-01 0. 


9.28235E-02 
9.30420E-02 
9.32575E-02 
9.34703E-02 
9.36800E-02 
9.38869E-02 
9.40908E-02 
9.429 17E-02 
9.44898E-02 
9.46852E-02 
9.48775E-02 
9.50668E-02 
9.5253 IE-02 
9.54365E-02 
9.56173E-02 
9.57947E-02 
9.59696E-02 
9.61415E-02 
9.63104E-02 
9.64765E-02 
9.66399E-02 
9.68002E-02 
9.69575E-02 
9.71123E-02 
9.72641E-02 
9.74132E-02 
9.75595E-02 
9.77027E-02 
9.78436E-02 
9.798 13E-02 
9.81163E-02 
9.82486E-02 
9.83785E-02 
9.85053E-02 
9.86296E-02 
9.87509E-02 
9.88698E-02 
9.89859E-02 
9.90996E-02 
9.92106E-02 
9.93191E-02 
9.94248E-02 
9.95279E-02 
9.96284E-02 
9.97266E-02 
9.98223E-02 
9.99156E-02 
1.00006E-01 
1.00094E-01 
1.00180E-01 
1.00264E-01 
1.00344E-01 
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0.584338 0. 
0.587944 0. 
0.591553 0. 
0.595159 0. 
0.598767 0. 
0.602373 0. 
0.605979 0. 
0.609588 0. 
0.613194 0. 
0.616802 0. 
0.620408 0. 
0.624014 0. 
0.627623 0. 
0.631229 0. 
0.634837 0. 
0.638443 0. 
0.642049 0. 
0.645658 0. 
0.649264 0. 
0.652872 0. 
0.656478 0. 
0.660084 0. 
0.663693 0. 
0.667299 0. 
0.670907 0. 
0.674513 0. 
0.678122 0. 
0.681728 0. 
0.685334 0. 
0.688942 0. 
0.692548 0. 
0.696157 0. 
0.699763 0. 
0.703369 0. 
0.706977 0. 
0.710583 0. 
0.714192 0. 
0.717798 0. 
0.721404 0. 
0.725013 0. 


1.00423E-01 0. 
1.00499E-01 0. 
1.00574E-01 0. 
1.00645E-01 0. 
1.00714E-01 0. 
1.00782E-01 0. 
1.00846E-01 0. 
1.00909E-01 0. 
1.00970E-01 0. 
1.01028E-01 0. 
1.01085E-01 0. 
1 .01 139E-01 0. 
1.01 191E-01 0. 
1.01241E-01 0. 
1.01289E-01 0. 
1.01335E-01 0. 
1.01379E-01 0. 
1 .0 1 42 1 E-0 1 0. 
1.01462E-01 0. 
1.01500E-01 0. 
1.01536E-01 0. 
1.0157 IE-01 0. 
1.01604E-01 0. 
1.01636E-01 0. 
1.01665E-01 0. 
1.01693E-01 0. 
1.01719E-01 0. 
1.01744E-01 0. 
1.01767E-01 0. 
1.01789E-01 0. 
1.01808E-01 0. 
1.01827E-01 0. 
1.01844E-01 0. 
1.01861E-01 0. 
1.01875E-01 0. 
1.01889E-01 0. 
1.01901E-01 0. 
1.01912E-01 0. 
1.01922E-01 0. 
1.01932E-01 0. 


1.00423E-01 
1.00499E-01 
1.00574E-01 
1.00645E-01 
1.00714E-01 
1.00782E-01 
1.00846E-01 
1.00909E-01 
1.00970E-01 
1.01028E-01 
1.01085E-01 
1.01139E-01 
1.01 191E-01 
1.01241E-01 
1.01289E-01 
1.01335E-01 
1.01379E-01 
1.01421E-01 
1.01462E-01 
1.01500E-01 
1.01536E-01 
1.01571E-01 
1.01604E-01 
1.01636E-01 
1.01665E-01 
1.01693E-01 
1.01719E-01 
1.01744E-01 
1.01767E-01 
1.01789E-01 
1.01808E-01 
1.01827E-01 
1.01844E-01 
1.01861E-01 
1.01875E-01 
1.01889E-01 
1.01901E-01 
1.01912E-01 
1.01922E-01 
1.01932E-01 


The following is the exact coordinates for the Long Mach 2.4 nozzle. The first 
line is minimum and maximum x values. Then the values are x, minimum y, maximum 
y, minimum z, and maximum z for one quadrant of the nozzle. Therefore, at any x 
location, the contour in the centerplane is found by taking the minimum z value and the 
maximum y value. 


-0.174889 1.29548 

-0.174889 0. 8.34442E-02 0. 8.34442E-02 
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-0.168412 0. 8.2195 IE-02 0. 8.21951E-02 
-0.161934 0. 8.09820E-02 0. 8.09820E-02 
-0.155457 0. 7.981 19E-02 0. 7.981 19E-02 
-0.148980 0. 7.86859E-02 0. 7.86859E-02 
-0.142502 0. 7.76042E-02 0. 7.76042E-02 
-0.136025 0. 7.65679E-02 0. 7.65679E-02 
-0.129548 0. 7.55770E-02 0. 7.55770E-02 
-0.123070 0. 7.463 15E-02 0. 7.463 15E-02 
-0.116593 0. 7.373 17E-02 0. 7.37317E-02 
-0.110115 0. 7.28783E-02 0. 7.28783E-02 
-1.03638E-01 0. 7.20709E-02 0. 7.20709E-02 
-9.71608E-02 0. 7.13101E-02 0. 7.13101E-02 
-9.06832E-02 0. 7.05956E-02 0. 7.05956E-02 
-8.42059E-02 0. 6.99278E-02 0. 6.99278E-02 
-7.77286E-02 0. 6.93073E-02 0. 6.93073E-02 
-7.12513E-02 0. 6.87336E-02 0. 6.87336E-02 
-6.47737E-02 0. 6.82069E-02 0. 6.82069E-02 
-5.82964E-02 0. 6.77278E-02 0. 6.77278E-02 
-5.18190E-02 0. 6.72962E-02 0. 6.72962E-02 
-4.53417E-02 0. 6.69124E-02 0. 6.69124E-02 
-3.88641E-02 0. 6.65765E-02 0. 6.65765E-02 
-3.23868E-02 0. 6.62885E-02 0. 6.62885E-02 
-2.59094E-02 0. 6.60486E-02 0. 6.60486E-02 
-1.94321E-02 0. 6.58572E-02 0. 6.58572E-02 
-1.29548E-02 0. 6.57146E-02 0. 6.57146&02 
-6.47740E-03 0. 6.56195E-02 0. 6.56195E-02 
0. 0. 6.55756E-02 0. 6.55756E-02 
6.47740E-03 0. 6.55856E-02 0. 6.55856E-02 
1.29548E-02 0. 6.56329E-02 0. 6.56329E-02 
1.94322E-02 0. 6.572 10E-02 0. 6.572 10E-02 
2.59095E-02 0. 6.58420E-02 0. 6.58420E-02 
3.23868E-02 0. 6.59886E-02 0. 6.59886E-02 
3.88644E-02 0. 6.61553E-02 0. 6.61553E-02 
4.534 17E-02 0. 6.63382E-02 0. 6.63382E-02 
5.18190E-02 0. 6.65333E-02 0. 6.65333E-02 
5.82964E-02 0. 6.67384E-02 0. 6.67384E-02 
6.47740E-02 0. 6.69509E-02 0. 6.69509E-02 
7.12513E-02 0. 6.71691E-02 0. 6.71691E-02 
7.77286E-02 0. 6.73910E-02 0. 6.73910E-02 
8.42059E-02 0. 6.76156E-02 0. 6.76156E-02 
9.06835E-02 0. 6.78418E-02 0. 6.78418E-02 
9.71608E-02 0. 6.80692E-02 0. 6.80692E-02 
1.03638E-01 0. 6.8297 IE-02 0. 6.8297 IE-02 
0.110115 0. 6.85257E-02 0. 6.85257E-02 
0.116593 0. 6.87547E-02 0. 6.87547E-02 
0.123070 0. 6.89832E-02 0. 6.89832E-02 
0.129548 0. 6.921 16E-02 0. 6.921 16E-02 
0.136025 0. 6.94395E-02 0. 6.94395E-02 
0.142502 0. 6.96675E-02 0. 6.96675E-02 
0.148980 0. 6.98958E-02 0. 6.98958E-02 
0.155457 0. 7.01238E-02 0. 7.01238E-02 
0.161934 0. 7.03518E-02 0. 7.03518E-02 
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0.168412 0. 7.05795E-02 0. 7.05795E-02 
0.174889 0. 7.08075E-02 0. 7.08075E-02 
0.181367 0. 7.10355E-02 0. 7.10355E-02 
0.187844 0. 7.12635E-02 0. 7.12635E-02 
0.194321 0. 7.1491 IE-02 0. 7.1491 IE-02 
0.200799 0. 7.17191E-02 0. 7.17191E-02 
0.207276 0. 7.19471E-02 0. 7.19471E-02 
0.213754 0. 7.21748E-02 0. 7.21748E-02 
0.220231 0. 7.24028E-02 0. 7.24028E-02 
0.226708 0. 7.26308E-02 0. 7.26308E-02 
0.233186 0. 7.28585E-02 0. 7.28585E-02 
0.239663 0. 7.30865E-02 0. 7.30865E-02 
0.246140 0. 7.33145E-02 0. 7.33145E-02 
0.252618 0. 7.35425E-02 0. 7.35425E-02 
0.259095 0. 7.37701E-02 0. 7.37701E-02 
0.265573 0. 7.3998 IE-02 0. 7.3998 IE-02 
0.272050 0. 7.42261E-02 0. 7.42261E-02 
0.278527 0. 7.44541E-02 0. 7.44541E-02 
0.285005 0. 7.4682 IE-02 0. 7.4682 IE-02 
0.291482 0. 7.49101E-02 0. 7.49101E-02 
0.297959 0. 7.51381E-02 0. 7.51381E-02 
0.304437 0. 7.53661E-02 0. 7.53661E-02 
0.310914 0. 7.55941E-02 0. 7.55941E-02 
0.317391 0. 7.58220E-02 0. 7.58220E-02 
0.323868 0. 7.60503E-02 0. 7.60503E-02 
0.330345 0. 7.62783E-02 0. 7.62783E-02 
0.336822 0. 7.65063E-02 0. 7.65063E-02 
0.343302 0. 7.67346E-02 0. 7.67346E-02 
0.349779 0. 7.69626E-02 0. 7.69626E-02 
0.356256 0. 7.71909E-02 0. 7.71909E-02 
0.362733 0. 7.74192E-02 0. 7.74192E-02 
0.369210 0. 7.76472E-02 0. 7.76472E-02 
0.375687 0. 7.78755E-02 0. 7.78755E-02 
0.382164 0. 7.81038E-02 0. 7.81038E-02 
0.388641 0. 7.83324E-02 0. 7.83324E-02 
0.395121 0. 7.85613E-02 0. 7.85613E-02 
0.401598 0. 7.87905E-02 0. 7.87905E-02 
0.408075 0. 7.90197E-02 0. 7.90197E-02 
0.414552 0. 7.92486E-02 0. 7.92486E-02 
0.421029 0. 7.94772E-02 0. 7.94772E-02 
0.427506 0. 7.97055E-02 0. 7.97055E-02 
0.433983 0. 7.99341E-02 0. 7.99341E-02 
0.440460 0. 8.01624E-02 0. 8.01624E-02 
0.446940 0. 8.03910E-02 0. 8.03910E-02 
0.453417 0. 8.06193E-02 0. 8.06193E-02 
0.459894 0. 8.08476E-02 0. 8.08476E-02 
0.466371 0. 8.10762E-02 0. 8.10762E-02 
0.472848 0. 8.13045E-02 0. 8.13045E-02 
0.479325 0. 8.15328E-02 0. 8.15328E-02 
0.485802 0. 8.17611E-02 0. 8.17611E-02 
0.492279 0. 8.19894E-02 0. 8.19894E-02 
0.498760 0. 8.22180E-02 0. 8.22180E-02 


0.505237 0. 
0.511714 0. 
0.518190 0. 
0.524668 0. 
0.531144 0. 
0.537621 0. 
0.544098 0. 
0.550579 0. 
0.557056 0. 
0.563533 0. 
0.570010 0. 
0.576487 0. 
0.582964 0. 
0.589441 0. 
0.595918 0. 
0.602398 0. 
0.608875 0. 
0.615352 0. 
0.621829 0. 
0.628306 0. 
0.634783 0. 
0.641260 0. 
0.647737 0. 
0.654217 0. 
0.660694 0. 
0.667171 0. 
0.673648 0. 
0.680125 0. 
0.686602 0. 
0.693079 0. 
0.699556 0. 
0.706036 0. 
0.712513 0. 
0.718990 0. 
0.725467 0. 
0.731944 0. 
0.738421 0. 
0.744898 0. 
0.751375 0. 
0.757855 0. 
0.764332 0. 
0.770809 0. 
0.777286 0. 
0.783763 0. 
0.790240 0. 
0.796717 0. 
0.803194 0. 
0.809674 0. 
0.816151 0. 
0.822628 0. 
0.829105 0. 
0.835582 0. 


8.24463E-02 0. 
8.26746E-02 0. 
8.29029E-02 0. 
8.31312E-02 0. 
8.33598E-02 0. 
8.35880E-02 0. 
8.38163E-02 0. 
8.40446E-02 0. 
8.42732E-02 0. 
8.45015E-02 0. 
8.47301E-02 0. 
8.49584E-02 0. 
8.51867E-02 0. 
8.54153E-02 0. 
8.56436E-02 0. 
8.58722E-02 0. 
8.61008E-02 0. 
8.6329 IE-02 0. 
8.65577E-02 0. 
8.67860E-02 0. 
8.70146E-02 0. 
8.72432E-02 0. 
8.747 18E-02 0. 
8.77004E-02 0. 
8.79287E-02 0. 
8.81573E-02 0. 
8.83859E-02 0. 
8.86145E-02 0. 
8.8843 IE-02 0. 
8.90720E-02 0. 
8.93006E-02 0. 
8.95292E-02 0. 
8.97578E-02 0. 
8.99867E-02 0. 
9.02153E-02 0. 
9.04439E-02 0. 
9.06728E-02 0. 
9.090 14E-02 0. 
9.11303E-02 0. 
9.13589E-02 0. 
9.15878E-02 0. 
9.18167E-02 0. 
9.20453E-02 0. 
9.22742E-02 0. 
9.25028E-02 0. 
9.273 17E-02 0. 
9.29600E-02 0. 
9.31883E-02 0. 
9.34163E-02 0. 
9.36437E-02 0. 
9.38705E-02 0. 
9.40963E-02 0. 


8.24463E-02 
8.26746E-02 
8.29029E-02 
8.313 12E-02 
8.33598E-02 
8.35880E-02 
8.381 63 E-02 
8.40446E-02 
8.42732E-02 
8.45015E-02 
8.47301E-02 
8.49584E-02 
8.51867E-02 
8.54153E-02 
8.56436E-02 
8.58722E-02 
8.61008E-02 
8.63291E-02 
8.65577E-02 
8.67860E-02 
8.70146E-02 
8.72432E-02 
8.747 18E-02 
8.77004E-02 
8.79287E-02 
8.81573E-02 
8.83859E-02 
8.86145E-02 
8.8843 IE-02 
8.90720E-02 
8.93006E-02 
8.95292E-02 
8.97578E-02 
8.99867E-02 
9.02153E-02 
9.04439E-02 
9.06728E-02 
9.09014E-02 
9.11303E-02 
9. 13589E-02 
9. 15878E-02 
9.18167E-02 
9.20453E-02 
9.22742E-02 
9.25028E-02 
9.273 17E-02 
9.29600E-02 
9.31883E-02 
9.34163E-02 
9.36437E-02 
9.38705E-02 
9.40963E-02 


0.842059 0 
0.848536 0 
0.855013 0 
0.861493 0 
0.867970 0 
0.874447 0 
0.880924 0 
0.887401 0 
0.893878 0 
0.900355 0 
0.906832 0 
0.913312 0 
0.919789 0 
0.926266 0 
0.932743 0 
0.939220 0 
0.945697 0 
0.952174 0 
0.958651 0 
0.965131 0 
0.971608 0 
0.978085 0 
0.984562 0 
0.991039 0 
0.997516 0 
1.00399 0. 
1.01047 0. 
1.01695 0. 
1.02343 0. 
1.02990 0. 
1.03638 0. 
1.04286 0. 
1.04934 0. 
1.05581 0. 
1.06229 0. 
1.06877 0. 
1.07525 0. 
1.08172 0. 
1.08820 0. 
1.09468 0. 
1.10115 0. 
1.10763 0. 
1.11411 0. 
1.12059 0. 
1.12707 0. 
1.13354 0. 
1.14002 0. 
1.14650 0. 
1.15297 0. 
1.15945 0. 
1.16593 0. 
1.17241 0. 


9.432 13E-02 0 
9.45453E-02 0 
9.47678E-02 0 
9.49882E-02 0 
9.52067E-02 0 
9.54228E-02 0 
9.56365E-02 0 
9.58474E-02 0 
9.60553E-02 0 
9.62604E-02 0 
9.64622E-02 0 
9.66606E-02 0 
9.68560E-02 0 
9.70480E-02 0 
9.72364E-02 0 
9.742 14E-02 0 
9.7603 IE-02 0 
9.77807E-02 0 
9.7955 IE-02 0 
9.81255E-02 0 
9.82922E-02 0 
9.84553E-02 0 
9.86147E-02 0 
9.87701E-02 0 
9.892 19E-02 0 
9.90698E-02 0. 
9.92139E-02 0. 
9.93538E-02 0. 
9.94904E-02 0. 
9.96230E-02 0. 
9.975 16E-02 0. 
9.98766E-02 0. 
9.99973E-02 0. 
1.001 14E-01 0. 
1.00228E-01 0. 
1.00337E-01 0. 
1.00443E-01 0. 
1.00545E-01 0. 
1.00643E-01 0. 
1.00738E-01 0. 
1.00829E-01 0. 
1.00916E-01 0. 
1.01000E-01 0. 
1.01080E-01 0. 
1.01156E-01 0. 
1.01229E-01 0. 
1.01298E-01 0. 
1.01365E-01 0. 
1.01427E-01 0. 
1.01486E-01 0. 
1.01542E-01 0. 
1.01595E-01 0. 


9.432 13E-02 

9.45453E-02 

9.47678E-02 

9.49882E-02 

9.52067E-02 

9.54228E-02 

9.56365E-02 

9.58474E-02 

9.60553E-02 

9.62604E-02 

9.64622E-02 

9.66606E-02 

9.68560E-02 

9.70480E-02 

9.72364E-02 

9.742 14E-02 

9.7603 IE-02 

9.77807E-02 

9.7955 IE-02 

9.81255E-02 

9.82922E-02 

9.84553E-02 

9.86147E-02 

9.8770 IE-02 

9.892 19E-02 

9.90698E-02 

9.92139E-02 

9.93538E-02 

9.94904E-02 

9.96230E-02 

9.975 16E-02 

9.98766E-02 

9.99973E-02 

1.001 14E-01 

1.00228E-01 

1.00337E-01 

1.00443E-01 

1.00545E-01 

1.00643E-01 

1.00738E-01 

1.00829E-01 

1.00916E-01 

1.01000E-01 

1.01080E-01 

1.01 156E-01 

1.01229E-01 

1.01298E-01 

1.01365E-01 

1.01427E-01 

1.01486E-01 

1.01542E-01 

1.01595E-01 
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1.17888 0. 1.01644E-01 0. 1.01644E-01 

1.18536 0. 1.01691E-01 0. 1.01691E-01 

1.19184 0. 1.01734E-01 0. 1.01734E-01 

1.19832 0. 1.01775E-01 0. 1.01775E-01 

1.20479 0. 1.01812E-01 0. 1.01812E-01 

1.21127 0. 1.01846E-01 0. 1.01846E-01 

1.21775 0. 1.01878E-01 0. 1.01878E-01 

1.22422 0. 1.01907E-01 0. 1.01907E-01 

1.23070 0. 1.01934E-01 0. 1.01934E-01 

1.23718 0. 1.01958E-01 0. 1.01958E-01 

1.24366 0. 1.01979E-01 0. 1.01979E-01 

1.25013 0. 1.01999E-01 0. 1.01999E-01 

1.25661 0. 1.02016E-01 0. 1.02016E-01 

1.26309 0. 1 .0203 IE-01 0. 1.02031E-01 

1.26957 0. 1.02044E-01 0. 1.02044E-01 

1.27604 0. 1.02055E-01 0. 1.02055E-01 

1.28252 0. 1.02065E-01 0. 1.02065E-01 

1.28900 0. 1.02073E-01 0. 1.02073E-01 

1.29548 0. 1.02081E-01 0. 1.02081E-01 


The program listed below is used for the generation of the grid for the Navier- 
Stokes calculation. It is called 3dgrid.f. Output from this file is binary data and is not 
shown here. 


* This program generates a three-dimensional nozzle grid for a N-S 

* calculation using an exponential stretching function adapted from 

* Dr. John Korte of NASA Langley. 

* 

* by: 

* 

* Tim Alcenius 

* Initially developed during LARSS Program participation 

* August 10, 1993 

* 

* Define Variables. 

* 

real xbeg, xend, delx, ymin(501), ymax(501), zmin(501), zmax(501) 
real sratio, beta, sexp, sfactc, xbar, ex, sratc, alp, sexc 
real xloc(501), x, y, z, sdec 
real r(305,97,97,3) 

integer choice, choice2, imax(l), jmax(l), kmax(l) 
integer i, j, k, 1, m 

* 

* Open input and output files and ask for stretching factor and 

* output format. 


* ^ 
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open(l, file = 'par am') 
open(2, file = 'geom', form = 'unformatted') 
print *, What stretching factor near the wall?' 
read *, beta 

5 print *, 'Do you want Tecplot data?' 
print *, ’1) yes' 
print *, ’0) no' 
read *, choice 

if(choice .NE. 1 .AND. choice .NE. 0) then 
goto 5 
endif 


* Open files for Tecplot format. 

* 

if (choice .EQ. 1) then 
open(3, file = 'grid') 
open(4, file = 'gridy') 
open(8, file = 'gridz') 
endif 

* 

* Set maximum array dimensions for all directions and set ymax 

* variable for interpolation routine. 

* 

imax(l) = 201 
jmax(l) = 65 
kmax(l) = 65 
do4i = 1,501 

ymax(i) = 0.000000000000000000000000000 
continue 

* Set all stretching constants. 

* 

xbar = . 1 
alp = 0 
s facte = 1 .25 
ex = 4. 
m = 1 

sratc = (sfactc +l.)/(sfactc - 1.) 
sratio = (beta+l.)/(beta-l.) 

* 

* Write initial Tecplot data and inital data in binary format 

* for N-S calculation. 

* 

if(chioce .E Q. 1) then 

write(3,*) T1 TlE=" S treamline Plot For Supersonic Nozzle’" 
write(3,*) ’VARIABLES="X","Y","Z’" 
endif 

write(2)(imax(l) jmax(l),kmax(l),l= 1,1) 

* 

* Read in data from file containing nozzle coordinates. 
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read(l,*) xbeg, xend 

delx = (xend-xbeg)/FLOAT(imax( 1 )- 1 ) 

x = xbeg 

do 10 i= 1,228 

read(l,*) xloc(i), ymin(i), ymax(i), zmin(i), zmax(i) 

10 continue 

do 20 i = 1, imax(l) 

* 

* Call interpolation routine for placing equally spaced points 

* in the x-direction independant of the spacing of data in the 

* input file. 

* 

* call interpolate(x,xloc,ymin,ymax,zmin,zmax,ymn,ymx,zmn,zmx,i) 

* Stretching functions in y and z directions. 

* 

do 30 j = 1, jmax(l) 
sexpy = (real(j-l)/real(jmax(l)-l) - alp ) 

+ / (1.0 -alp) 

sdeny = 1. + sratio**sexpy 
sexcy = real(jmax(l)-j)/real(jmax(l)-l) 
sdecy = 1. + sratc**sexcy 
y = ymn 

+ xbar * sexcy* *ex 

* ( ( ymx - ymn ) 

* ( 1.0 - sfactc + 2.0*sfactc/sdecy ) ) 

+ ( ( 1. - xbar ) + xbar * ( 1. - sexcy**ex ) ) 

* (ymx - ymn ) 

* (2.0* alp/ (2.0* alp + 1) 

+ beta / (2.0 * alp + 1) * (1.0 - 2.0/sdeny) ) 
do 40 k = 1, kmax(l) 
sexpz = (real(k-l)A-eal(kmax(l)-l) - alp ) 

+ / (1.0 -alp) 

sdenz = 1. + sratio**sexpz 
sexcz =real(kmax(l)-k)/real(kmax(l)-l) 
sdecz = 1. + sratc**sexcz 
z = zmn 

+ xbar * sexcz* *ex 

* ( ( zmx - zmn ) 

* ( 1.0 - sfactc + 2.0*sfactc/sdecz ) ) 

+ ( ( 1. - xbar ) + xbar * ( 1. - sexcz**ex ) ) 

* (zmx - zmn ) 

* (2.0* alp/ (2.0* alp + 1) 

+ beta / (2.0 * alp + 1) * ( 1 .0 - 2.0/sdenz) ) 

* 

* Write certain planes if Tecplot data is desired. 

* 

if(choice .EQ. 1) then 
if(m .EQ. i) then 
if(flagx .EQ. 0) then 

write(3,*) 'ZONE T="Wall Countour", I- jmax(l),’,J=', 

+ kmax(l),'J 7 =POINT 



137 


flagx = 1 
endif 

write(3,99) x, y, z 
endif 

if(y .EQ. ymx) then 
if(flagy EQ. 0) then 

write(4,*) 'ZONE T="Wall Countour",I=',kmax(l),' J=', 
+ imax(l),',F=POINT 

flagy = 1 
endif 

write(4,99) x, y, z 
endif 

if(z .LT. .0001 .AND. z .GT. 0) then 
if(flagz .EQ. 0) then 

write(8,*) 'ZONE T="Wall Countour",I=’jmax(l),',J=’, 
+ imax( l),',F=POINT 

flagz = 1 
enchf 

write(8,99) x, y, z 
endif 

if(i .GT. m) then 
m = m+40 
flagx = 0 
endif 
endif 

♦ 

* Put data in one array for use to write data in binary format. 

* 

r(i,j,k,l) = x 
r(i,j,k,2) = y 

r(i,j,k,3) = z 
40 continue 
30 continue 
x = x + delx 

20 continue 

* 

* Write binary file of data for N-S calculation. 

* 

write(2) ((( ((r(i,j ,kj)),i= 1 dmax( 1 )) ,j= 1 jmax( 1 )),k= 1 , 

+ kmax(l)),l=l,3) 

99 format(3(lx,el4.7)) 

end 


* Interpolation subroutine. 

* 

subroutine interpolate^, xloc,ymin,ymax,zmin,zmax,ymn,ymx,zmn, 
+ zmx,i) 
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* Define variables. 

* 

real xloc(501), ymin(501), ymax(501),zmin(501), zmax(501) 
integer ije 

* 

* Set dummy variable so counter in program is unchanged. 

* 

k = i 

* 

* Make sure that counter does not exit the boundaries of the 

* nozzle parameter file. 

* 

100 if(ymax(k) .EQ. 0) then 
k = k-1 
goto 100 
endif 

* 

* If the x location corresponds with the parameter file then use 

* the y and z values for that location. 

* 

if(ABS(xloc(k)-x) .LT. .00001) then 
200 ymn = ymin(k) 
ymx = ymax(k) 
zmn = zmin(k) 
zmx = zmax(k) 
return 
endif 

* 

* Make sure that counter does not exit the boundaries of the 

* nozzle parameter file. 

* 

if( k .LT. 1) then 
k = k+1 
goto 200 
endif 

* 

* Find boundaries on x and determine y and z by using a first 

* order interpolation routine. 

* 

if(xloc(k) .LT. x) then 
if(xloc(k+l) .GT. x) then 

ymn = ((ymin(k+l)-ymin(k))*(x-xloc(k)))/(xloc(k+l)-xloc(k))+ 
+ ymin(k) 

ymx = ((ymax(k+l)-ymax(k))*(x-xloc(k)))/(xloc(k+l)-xloc(k))+ 
+ ymax(k) 

zmn = ((zmin(k+l)-zmin(k))*(x-xloc(k)))/(xloc(k+l)-xloc(k))+ 
+ zmin(k) 

zmx = ((zmax(k+l)-zmax(k))*(x-xloc(k)))/(xloc(k+l)-xloc(k))+ 
+ zmax(k) 
return 
else 
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k = k+1 
goto 100 
endif 
else 

* 

* Find boundaries on x and determine y and z by using a first 

* order interpolation routine. 

* 

if(xloc(k-l) .LT. x) then 

ymn = (ymin(k)-ymin(k-l))*(x-xloc(k-l))/(xloc(k)-xloc(k-l))+ 
+ ymin(k-l) 

ymx = (ymax(k)-ymax(k-l))*(x-xloc(k-l))/(xloc(k)-xloc(k-l))+ 
+ ymax(k-l) 

zmn = (zmin(k)-zmin(k-l))*(x-xloc(k-l))/(xloc(k)-xloc(k-l))+ 
+ zmin(k-l) 

zmx = (zmax(k)-zmax(k-l))*(x-xloc(k-l))/(xloc(k)-xloc(k-l))+ 
+ zmax(k-l) 
return 
else 
k = k-1 
goto 100 
endif 
endif 
return 
end 


A ppendix B Post-Processing Code 

The program below was developed to calculate the crossflow from the results of 
the numerical simulation. 

* This program has been made to calculate various flowfield 

* parameters after a three-dimensional N-S nozzle calculation. 

* 

* by: 

* 

* Tim Alcenius 

* Purdue University 

* August 13, 1994 

* 

* Identify variables 

* 

real mach(305,97,97), p(305,97,97), wmaxpct(305,97) 

real r(98,98,306,5), x(305,97,97,3) 

real uold, void, wold, pOold, diffp(305), t2, tl 

real unew, vnew, wnew, pint, machint, pOnew, ymax, capu, pO 

real si, s2, s3, sdiv, ndistnew, prof(97,6), capuml 

real magV, shpfct(305,97), nl, n2, n3, ql, q2, q3 



real qdiv, delst(305,97), theta(305,97), c5, c2, c3, c4 
real grtheta, beta, xxsi, yxsi, zxsi, xeta, yeta, zeta 
real rhoe, nue, crss, pOl, p02, Hint, tdellO, gamma, pi 
real Ve, magVo, wmax, magucf, dellO, H, rcf(305,97) 
real tnew, told, xnew, xold, znew, zold, rcfnew(305,97) 
real A, L, cst, mache, cstad, tstote, ndist, ndisttot 
real Cl(305,97), rcfd2int, rcfdel2(305,97), tstadote 
real rhoold, rhonew, nminus, magucfo 
integer imax(l), jmax(l), kmax(l), i, j, k, m 
integer iref, kedge, jedge, flag 
pi = ACOS(-l.) 

Open needed files and read data from LARCK output files written 
in binary format. 

open(l, file = 'geom', form = 'unformatted') 

open(2, file = 'plt3d.q', form = 'unformatted') 

open(7, file = 'diffp') 

open(8, file = 'streammag') 

open(9, file = 'bledge') 

open(10, file = 'deist') 

open(l 1, file = 'theta') 

open(14, file = 'crossflow') 

open(15, file = 'comer') 

open(16, file = 'crossflowmag') 

read( 1 )(imax(m), jmax(m), kmax(m), m=l,l) 

read( l)((((x(i j.kjn), i=l,imax(l)), j=l,jmax(l)), k=l,kmax(l)), 

+ m=l,3) 

read(2)(imax(m), jmax(m), kmax(m), m=l,l) 
read(2) c5, c2, c3, c4 

read(2)((((r(j,k,i,m), i=l,imax(l)), j=l jmax(l)), k=l,kmax(l)), 

+ m=l,5) 

imax(l) = imax(l) - 1 
jmax(l) = jmax(l) - 1 
kmax(l) = kmax(l) - 1 

Set-up Tecplot format in top of output files. 

write(7,*)''nTLE= "Delta P’" 

write(7,*)’VARIABLES="x", "delta P"’ 

write(7,*)'ZONE T="planes",I=',imax( 1 )- l,',F=POINT 

write(8,*)TlTLE="Streamwise'" 

write(8,*) r VARIABLES="y","V"' 

write(9,*)TirLE="Streamwise"' 

write(9,*)'VARIABLES="y","edge"' 

write(10,*)'TITLE="delst'" 

write(10,*)'VARIABLES="x", "deist'" 

write( 10,*)'ZONE T="planes",I=’,imax(l)-l,',F=POINT 

write( 1 1 *)"nTLE="theta'" 

write(ll,*)'VARIABLES="x", "theta"' 

write(l 1 ,*)'ZONE T="planes",I='4max(l)-l,'^=POINT 

write( 1 4, *)'TITLE="Crossflow Re’" 



write(14,*)'VARIABLES="x","y"»"Conventional","Reed&Haynes", 

+ "Wmax/Ue”,"R","Amal","Arnal Cr',"Shp"' 

write(14,*)'Z0NE T= "planes", I=',jmax(l)-2,' ,J=',imax(l)- 

+ l,'J T =POrNT 

write(15 *) , TrrLE="Corner"' 

write(15,*)'VARIABLES="x","y","u","v","mach’" 

write( 1 5,*) , Z0N E T="planes",I=',jmax( 1 ),',J=',imax( 1 )- 1 ,',F=POINT 

write(16,*)'TlTLE="Crossflow'" 

write(16,*)'VARIABLES="y","w"’ 


* 

* Calculate variables u, v, w from conserved variables rho*u, 

* rho*v, and rho*w. Also pressure, Mach number, 

* and total pressure. 

* 

do 10 i= l,imax(l) 
do 20 j = ljmax(l) 
do 30 k = l,kmax(l) 
do 40 m = 2,4 

if(m .EQ. 2) r(j,k,i,2) = r(j,k,i,2)/r(jjc,i,l) 
if(m .EQ. 3) r(j,k4,3) = r(j,k,i,3)/r(j,k,i,l) 
if(m .EQ. 4) r(j,k,i,4) = r(j,k,i,4)/r(j,k,i,l) 

40 continue 

p(i,j,k) = .4*(r(j,k,i,5)-.5*r(j,k,i,l)*(r(jjc,i,2)**2+ 

+ r(j,k,i,3)**2+r(j,k,i,4)**2)) 

mach(i,j,k) = SQRT((r(j,k,i,2)**2+r(j,k,i,3)**2+ 

+ r(j,k4,4)**2)/(1.4*p(i,j,k)/r(j,k4,l))) 

30 continue 

wmaxpct(i.j) = -.0005 
20 continue 
10 continue 

pO = p(2, 1 , 1 )*( 1 .+.2*mach(2, 1 , 1 )**2)**(3.5) 

* 

* Call subroutine to outer limit of marching in any 

* constant i-plane. 

* 

do 150 i = 2,imax(l) 

call findjedge(p0, jedge, mach, p, i, jmax) 
diffp(i) = (p(i jmax( 1 ),kmax( 1 ))-p(i, 1 ,kmax( 1 )))^>0 
write(7,*) x(i,l,l,l),diffp(i) 

♦ 

* Calculate the inward normal to the wall 

* 

do 200 j = 1, jedge 
iref = i 

xxsi = x(i+l j,kmax(l),l)-x(ij,kmax(l),l) 
yxsi = x(i+l,j,kmax(l),2)-x(i,j,kmax(l),2) 
zxsi = x(i+l,j,kmax(l),3)-x(i,j,kmax(l),3) 
xeta = x(i,j+l,kmax(l),l)-x(i,j,kmax(l),l) 
yeta = x(i,j+ldcmax(l),2)-x(i,j,kmax(l),2) 
zeta = x(i,j+l,kmax(l),3)-x(i,j,kmax(l),3) 
nl = (yeta*zxsi-yxsi*zeta)/(SQRT((yeta*zxsi-yxsi*zeta)**2+ 
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+ (xeta*zxsi-xxsi*zeta)**2+(xxsi*yeta-yxsi*xeta)**2)) 
n2 = -nl*((xeta*zxsi-xxsi*zeta)/(yeta*zxsi-yxsi*zeta)) 
n3 = (-xxsi*nl-yxsi*n2)/zxsi 

* 

* Set inital variables and reset all integral counters. 

* 

H = 0.0 
Hint = 0.0 
wmax = 0.0000 
dellO = 0. 
delst(i,j) = 0. 
theta(ij) = 0. 
rcfd2int = 0. 

told = p(i,j,kmax(l))/(287.*r(j,kmax(l),i,l)) 

xold = x(i,j,kmax(l),l) 
zold = x(ij,kmax(l),3) 
uold = r(j,kmax(l),i,2) 
void = r(j,kmax(l),i,3) 
wold - r(j,kmax(l),i,4) 
rhoold = r(j,kmax(l),i,l) 

p0old= p(i,j,kmax(l))*(l.+.2*mach(i,j,kmax(l))**2)**(3.5) 
ndisttot = 0.0 

* 

* Check to make sure normal vector doesn't cross a line outside 

* of the domain. 

* 

if(nl .LT. 0) then 
iref = i - 1 
else 

iref = i + 1 
endif 

if (iref .LE. 0 .OR. iref .GT. imax(l)) then 
wmax = 0.0 
goto 1010 
endif 
flag = 0 

* 

* Find each location that the normal vector crosses a z = constant 

* line to find wmax, . 1 *wmax, and integrate to find H. 

do 220 k = kmax(l), 1, -1 
grtheta = ATAN(ABS(x(iref,j,k,3)-x(ij,k,3))/ 

+ ABS (x(iref,j Jc, 1 )-x(i j ,k, 1 ))) 

beta = ATAN(ABS(nl/n3)) 
gamma = pi/2. - beta + grtheta 
delx3 = (x(i,j,kmax(l),3)-x(ij,k,3))*SIN(beta)/ 

+ SEN(gamma) 

1000 diag = SQRT((x(irefj,k,3)-x(ij,k,3))**2+ 

+ (x(irefj,k,l)-x(i,j,k,l))**2) 

* 

* Make sure normal vector crosses grid lines inside the domain. 



if(delx3 .GT. diag) then 
if(nl .LT. 0) then 
iref = iref - 1 
else 

iref = iref + 1 
endif 

if(iref .LE. 0 .OR. iref .GT. imax(l)) then 
wmax = 0.0 
goto 1010 
endif 
goto 1000 
endif 

Calculate normal distance in grid cell for integration using 
trapezoidal method along line normal to the wall. Interpolate 
variables at each crossing and save temperature, distance, 
and velocity components. Finished when normal vector 
crosses 98% of Po on the centerline. 

znew = x(i j,k,3)+delx3*(x(iref j,k,3)-x(i,j,k,3))/diag 

xnew = x(ij,k,l)-i-delx3*(x(iref,j,k,l)-x(i,jjc,l))/diag 

t2 = p(iref,j,k)/(287 .*r(j,k,iref, 1 )) 

tl =p(ij,k)/(287.*r(j,k,i,l)) 

tnew = (tl + (t2-tl)*delx3/diag) 

ndist = SQRT((xnew-xold)**2+(znew-zold)**2) 

unew = (r(j,k,i,2) + (r(j,k,iref,2)-r(j,k,i,2))* 

+ delx3/diag) 

vnew = (r(j,k,i,3) + (r(j,k,iref,3)-r(j,k,i,3))* 

+ delx3/diag) 

wnew = (r(j,k4,4) + (r(j,k,iref,4)-r(j,k,i,4))* 

+ delx3/diag) 

pint = (p(ij,k) + (p(irefj,k)-p(i,j,k))* 

+ delx3/diag) 

machint = (mach(ijjc) + (mach(iref,j,k)-mach(ij,k))* 

+ delx3/diag) 

rhonew = pint/(287.*tnew) 
p0new= pint*(l .+.2*machint**2)**(3.5) 
if(p0new .GT. .98*p0) then 
ndistnew = (.98*p0-p0old)*ndist/(p0new-p0old) 
ndisttot = ndisttot + ndistnew 
prof(k,l) = told + ndistnew*(tnew-told)/ndist 
prof(k,2) = ndisttot 

prof(k,3) = uold + ndistnew*(unew-uold)/ndist 
prof(k,4) = void + ndistnew*(vnew-vold)/ndist 
prof(k,5) = wold + ndistnew*(wnew-wold)/ndist 
prof(k,6) = rhoold + ndistnew*(rhonew-rhoold)/ndist 
mache = SQRT((prof(k,3)* !|t 2+prof(k,4)**2+prof(k^)**2)/ 
+ (1 .4*287*prof(k, 1))) 

rhoe = .98*p0/(287*prof(k,l)*(l.+.2*mache**2)**(3.5)) 
nue = (.17 16e-4*(prof(k, l)/273. 1 1 1)** 1 ,5*(383.444/ 

+ (prof(k,l)+110.333)))/rhoe 

kedge = k 
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goto 1020 
else 

ndisttot = ndisttot + ndist 
prof(k,l) = tnew 
prof(k,2) = ndisttot 
prof(k,3) = unew 
prof(k,4) = vnew 
prof(k,5) = wnew 
prof(k,6) = rhonew 
endif 

* 

* Set all old variables to current values for next step in 

* do-loop. 

* 

xold = xnew 
zold = znew 
told = tnew 
uold = unew 
void = vnew 
wold = wnew 
rhoold = rhonew 
pOold = pOnew 

220 continue 

* 

* Calculate crossflow direction unit normal. 

* 

1020 if((prof(kedge,3)*n3-prof(kedge,5)*nl) .NE. 0.) then 
q2 = 1 . 

q3 = (-prof(kedge,3)*n2+nl*prof(kedge,4))/(prof(kedge,3)* 
+ n3-prof(kedge,5)*nl) 

ql = -(prof(kedge,4)+q3*prof(kedge,5))/prof(kedge,3) 
qdiv = SQRT(ql**2+q2**2+q3**2) 
ql = ql/qdiv 
q2 = q2/qdiv 
q3 = q3/qdiv 
else 
ql =0. 
q2 = 1. 
q3 = 0. 
endif 

* 

* Calculate streamiwse unit vector. 

* 

si = 1. 

s3 = (ql *n2-q2*n 1 )/(n3 *q2-q3*n2) 

s2 = -(s3*q3+ql)/q2 

sdiv = SQRT(sl**2+s2**2+s3**2) 

si = sl/sdiv 

s2 = s2/sdiv 

s3 = s3/sdiv 

* 

* Find Wmax, del 10, and integrate for H, die*. Theta, and 



Rdel2 form Amal AGARD R-786. 

magucfo = 0. 
magVo = 0. 

Ve = SQRT(prof(kedge,3)**2+prof(kedge,4)**2+ 

+ prof(kedge,5)**2) 
do 260 k = kmax(l), kedge, -1 
magucf = prof(k,3)*ql + prof(k,4)*q2 + prof(k,5)*q3 
magV = prof(k,3)*sl + prof(k,4)*s2 + prof(k,5)*s3 
if(ABS (magucf) .GT. ABS(wmax) .AND. k .NE. kmax(l)) then 
wmax = magucf 
else 

if(k .EQ. kedge .AND. magucf*magucfo .LT. 0.) then 
magucf = -magucf 
endif 

if(ABS (magucf) .LE. . l*ABS(wmax) .AND. magucf* 

+ magucfo .GT. 0. .AND. ABS(magucfo) .GT. .1* 

+ ABS(wmax)) then 

if ( AB S (magucfo-magucf) .LT. .000001) then 
nminus =0. 
else 

nminus = (.l*wmax-magucf)*(prof(k,2)-prof(k+l,2))/ 

+ (magucfo-magucf) 

endif 

del 10 = prof(k,2) - nminus 

tdellO = prof(k+l,l) + (prof(k,l)-prof(k+l,l))* 

+ (dell0-prof(k+l ,2))/(prof(k,2)-prof(k+l,2)) 

H = Hint + .5*(dellO-prof(k+l,2))*(tdellO+ 

+ prof(k+l,l))/prof(kedge,l) 

endif 
endif 

Write crossflow profiles approximately half-way between 
comer and centerplane of the nozzle. 

if(j .EQ. 19) then 
if(flag .EQ. 0) then 

write(16,*)'ZONET="planes",I-,kmax(l)-kedge+l, 

+'T J =POENT 

write(8,*)'ZONE T="planes",I=',kmax(l)-kedge+l, 
+'JF=POINT 

write(9,*) x(ij,kedge,2),prof(kedge,2) 
flag = 1 
endif 

write (16,*) prof (k,2)/prof(kedge,2) .magucf 
write(8, * ) prof(k,2)/prof (kedge,2) jnag V 
endif 

Hint = Hint + .5*(prof(k-l,2)-prof(k,2))* 

+ (prof(k, l)+prof(k- 1 , 1 ))/prof(kedge, 1) 

if(k .NE. kmax(l)) then 

rcfd2int = rcfd2int + .5*(prof(k,2)-prof(k+l,2))* 

+ (ABS(magucf)+ABS(magucfo)) 



* 

* 


Calculate delta* and theta. 


capu = 1 . - (pr of(k,6)/pr of (kedge,6)) *(magV/Ve) 
capuml = l.-(prof(k+l,6)/prof(kedge,6))*(mag Vo/Ve) 
delst(i.j) = delst(i j) + .5 * (prof (k,2) -prof (k+ 

+ l,2))*(capuml+capu) 

capu = (l.-capu)*(l.-(magV/Ve)) 
capuml = (l.-capuml)*(l.-(magVo/Ve)) 
theta(ij) = theta(ij) + ,5*(prof(k,2)-prof(k+l,2))* 

+ (capuml +capu) 

endif 

magucfo = magucf 
magVo = magV 
260 continue 


* Calculate standar crossflow Reynolds number and crossflow 

* Reynolds number from Amal. 

* 


1010 rcf(ij) = ABS(wmax)*dellO/nue 
rcfdel2(ij) = rcfd2int/nue 
shpfct(ij) = delst(i,j)/theta(ij) 
if(shpfct(i,j) .GT. 2.3) then 

Cl(ij) = 300.*ATAN(.106/((shpfct(ij)-2.3)**2.05))/pi 

else 

if(shpfct(i,j) .EQ. 2.3) then 
Cl(i,j) = 150. 
else 

Cl(iJ)=0. 

endif 

endif 

* 

* Calculate crossflow Reynolds number according to eons 5-12 

* AIAA 93-3054 by Reed and Haynes. 


* 


A = SQRT(0.72)*0.2*mache**2 

tstote = 0.5 + 0.5*prof(kmax(l), l)/prof(kedge,l) + A/6 

tstadote = 0.5 + 0.5*(1.+A) + A/6. 

cst = (S QRT (ts tote)* ( 1+11 0.4/prof(kedge, 1 )))/(tstote+ 
+ 1 10.4/prof(kedge, 1)) 

cstad = (SQRT(tstadote)*(l+110.4/prof(kedge,l)))/ 

+ (tstadote+1 10.4/prof(kedge, 1)) 
if(H .EQ. 0.) then 


H = 1. 
else 


H = dell0/H 
endif 


L = (SQRT(cst/cstad)*(3. 279+1 .72 1 *( 1 ,+A)+0.664* A))/ 
+ (5.+2.385*A) 

rcfnew(ij) = H*L*rcf(ij) 


* Calculate Wmax/Ue in % and set 0. to slightly greater 



* than 0. for use in divsor later. 

* 

wmaxpct(i,j) = 100*wmax/SQRT(prof(kedge,3)**2+ 

+ prof(kedge,4)**2+prof(kedge,5)**2) 

if(ABS(wmaxpct(i,j)) .LT. 0.00005) then 
wmaxpct(i.j) = -.00005 
endif 

if(j .EQ. 1) then 
write(10,*) x(i,j,k,l), delst(ij) 
write(ll,*) x(i,j,k,l), theta(ij) 
endif 

200 continue 
* 

* Write files for Reynolds number contours and diagonal 

* plane variables. 

* 

do 230 m = 2, jmax( 1 )- 1 

write(14,*) x(i,m,kmax(l),l),x(i4n,kmax(l),2)/cf(i^n), 

+ rcfnew(i,m),wmaxpct(i,m),rcfnew(i,m)/ABS(wmaxpct(i,m)) > 

+ rcfdel2(i,m),Cl(i,m),shpfct(i,m) 

230 continue 

do 240 j = 1, jmax(l) 
do 250 k * 1, kmax(l) 
if(j .EQ. k) then 

write(15,*) x(i,j,k 1 l) ( SQRT(2.)*x(ij,k,2),r(j,k,i,2), 

+ SQRT(r(j,k,i,3)**2+r(j,k,i,4)**2),mach(i,j,k) 

endif 

250 continue 
240 continue 
150 continue 
end 


* This subroutine finds the boundary layer edge by determining 

* where the total pressure drops to 98% of its value on the 

* centerline of the nozzle. 

* 

* Define variables 

* 

subroutine findjedge(p0, jedge, mach, p, i, jmax) 
real locpO, pO, mach(305,97,97), p(305,97,97) 
integer flag, jedge, jmax(l), i, j 

* 

* Calculate local pO. If below 98% of the centerline value, 

* jedge is saved. 

* 

flag = 0 

do 20 j = 1, jmax(l) 

locpO = p(ij,l)*(l.+.2*mach(i,j,l)**2)**(3.5) 



if(locpO .LT. ,98*pO AND. flag .NE. 1) then 
jedge = j-1 
flag = 1 
endif 

20 continue 
return 
end 




